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ABSTRACT 

The  scattering  of  sound  by  the  nonlinear  interaction  of  two 
sound  beams  in  the  presence  of  turbulence  is  used  to  experimentally 
measure  the  turbulent  velocities  generated  by  a  submerged  water 
jet.  When  two  sound  beams  of  primary  frequencies  foi  andfo2 
intersect  in  a  region  of  turbulent  flow,  the  nonlinear  scattering 
generates  sum  and  difference  frequency  components  (fo+  =  foi  +102 
and  fo-  =  foi  -  fo2)  which  radiate  outside  the  interaction  region.  In 
the  absence  of  turbulence,  the  crossed  beams  do  not  produce 
radiated  sum  and  difference  frequencies.  In  this  experiment,  two 
transducers  emit  continuous  wave  focused  sound  beams  of 
frequencies  foi  =  2.0  MHz  andfo2  =  2.1  MHz,  respectively.  The  sound 
beams  are  arranged  so  that  their  focal  points  overlap  and  the  beam 
axes  are  mutually  perpendicular.  A  receiving  transducer,  located 
outside  the  interaction  region,  detects  the  scattered  sum  frequency 
(f0+  =  4.1  MHz). 

The  turbulent  velocities  in  the  small  interaction  volume  are 
determined  from  variations  in  the  shape  of  the  scattered  sum 
frequency's  intensity  spectrum  that  are  measured  at  each  scattering 
angle.  The  motion  of  the  turbulent  eddies  (which  are  responsible  for 
the  nonlinear  scattering)  generates  a  random  Doppler  shift  of  the 
sum  frequency  component  which  broadens  its  intensity  spectrum. 
Measurements  of  the  average  Doppler  shift,  rms  frequency, 
skewness,  and  kurtosis  of  the  time  averaged  spectra  are  used  to 
correlate  the  mean  and  turbulent  velocity  components  along  the 
radial  and  axial  directions  of  the  jet.  Scattering  results  measured  at 
80  angles  for  each  of  11  scan  positions  across  the  jet  are  used  to  map 
out  the  velocity  correlation  coefficients  of  the  turbulence.  Results  are 
then  compared  with  earlier  published  experiments  that  use 
conventional  hot  wire  probes  to  measure  a  similar  turbulent  jet  flow 
in  air. 
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HISTORICAL  INTRODUCTION 


This  underwater  experiment  studies  a  nonlinear  acoustic 
scattering  phenomenon  which  generates  sum  and  difference 
frequencies  as  a  result  of  the  scattering  of  sound  by  sound  in 
the  presence  of  turbulence.  It  is  well  known  that  when  two 
collinear  finite-amplitude  sound  beams  of  different  frequencies 
interact,  they  produce  two  new  frequencies  which  are  equal  to 
the  algebraic  sum  and  difference  of  the  original  two 
frequencies.  These  are  commonly  referred  to  as  "combination" 
frequencies.  When  these  same  two  sound  beams  (each 
considered  to  be  well  collimated)  interact  at  right  angles  to 
each  other,  the  combination  frequencies  are  again  created  but 
do  not  radiate  beyond  the  bounds  of  the  region  defined  by  the 
overlap  of  the  two  sound  beams1.  Fig.  1  graphically  describes 
the  overlap  region. 


beam  *1 


beam  *2 


FIG.  1.  Region  defined  by  overlapping  sound  beams 

In  order  to  produce  a  situation  where  the  sum  and  difference 
frequencies  can  radiate  outside  the  overlap  region  (of  the 


crossed  beams),  one  must  introduce  a  solid  object,  air  bubbles, 
inhomogeneous  medium,  or  a  turbulent  flow  in  the  interaction 
region10*11.  In  these  cases,  the  nonlinear  interaction  will  allow 
sum  and  difference  frequency  components  to  radiate. 

The  intent  of  this  paper  is  to  explore  the  nonlinear 
scattered  sum  frequency  component  when  turbulence  is 
present  in  the  interaction  region  of  the  crossed  sound  beams. § 
Experimental  measurements  of  the  characteristic  properties  of 
the  sum  frequency  component's  intensity  spectrum  will  be 
used  to  determine  many  statistical  quantities  associated  with 
the  turbulent  velocity  field. 

In  the  early  1980s,  Korman  and  Beyer4*5  developed  an 
experiment  to  study  the  nonlinear  scattering  of  crossed  sound 
beams  in  the  presence  of  turbulence.  Their  experimental 
apparatus  generated  mutually  perpendicular  acoustic  beams 
(each  with  half  power  beam  widths  of  about  1  degree)  from 
two  unfocused  circular  plane  array  transducers.  A  submerged 
circular  water  jet  created  the  turbulence  in  their  water  tank. 
With  this  apparatus,  they  measured  the  intensity  of  the 
scattered  sum  frequency  as  a  function  of  angle.  To  analyze  the 
data,  Korman  and  Beyer  modeled  the  turbulence  as  an  isotropic 
distribution  of  homogeneous  randomly  fluctuating  velocities 


§  Note  that  throughout  this  paper,  extensive  references  are 
made  to  the  sum  frequency  and  very  little  is  said  about  the 
difference  frequency  component.  This  is  not  an  oversight.  It  is 
known  that  the  relative  intensity  of  the  difference  frequency 
component  to  the  sum  frequency  component  is  very  small.  This 
makes  detection  of  the  difference  frequency  difficult.  Since  the 
wavelengths  of  the  sum  frequency  are  much  smaller  than  the 
difference  frequency,  one  can  perform  the  experiments  in  a 
small  laboratory  tank.  Further,  an  analysis  of  scattering  at  the 
difference  frequency  reveals  that  angular  scattering  is  very 
insensitive  to  the  characteristic  length  scales  of  the  turbulence 
in  contrast  to  scattering  at  the  sum  frequency. 
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superimposed  upon  a  mean  flow  velocity.  They  modeled  the 
turbulent  velocity  fluctuations  to  behave  with  Gaussian 
statistics.  With  this  model,  they  could  predict  the  average 
features  of  the  turbulent  flow  across  their  jet  from  the  shape  of 
the  scattered  acoustic  intensity  spectra.  Specifically,  the 
Doppler  shift,  cod,  of  the  sum  frequency  was  related  to  the  jet's 
mean  velocity,  V,  and  the  spectral  broadening  was  related  to 
the  rms  turbulent  flow  velocity,  o.  The  intensity  spectrum 
I+(f,  0,)  integrated  over  all  frequencies  in  a  band  near  the  sum 

frequency  component  I+(6.)  =  f  *+(*•  ^ 

was  related  to  the  turbulent  kinetic  energy  spectrum,  E(k), 
which  describes  the  wavenumber,  k,  of  the  turbulent  eddies. 
Here,  I+(0.)  is  the  total  intensity  at  the  scattering  angle  0.. 

In  1987,  Korman  worked  with  Rife  on  a  USNA  Trident 
Project8  to  measure  turbulent  velocities  with  focused  sound 
beams.  Focusing  the  beams  minimized  the  volume  of  the 
interaction  region  and  allowed  Korman  and  Rife  to  spatially 
resolve  the  velocity  information  derived  from  the  scattered 
spectra.  This  was  not  possible  in  Korman  and  Beyer’s  earlier 
work  because  the  volumes  defined  by  the  overlapping  beams 
were  too  large  to  resolve  any  details  of  the  jet. 

Korman  and  Rife6*7  performed  two  different  acoustic 
scattering  experiments  with  the  focused  ultrasonic  beams.  The 
first  involved  two  crossed  focused  beams  which  were 
translated  across  the  width  of  the  jet,  as  shown  in  Fig.  2. 


6 


FIG  2.  Geometry  of  translational  scans 

The  second  experiment  involved  the  "conventional"  scattering 
of  one  of  the  individual  primary  sound  beams  by  the 
turbulence.  Here,  the  scattered  primary  frequency’s  intensity 
spectrum  is  investigated.  This  experiment  was  identical  in 
geometry  to  the  first,  except  only  one  sound  beam  was 
translated  across  the  turbulent  jet.  In  both  experiments,  the 
relative  angles  for  both  of  the  transmitter  beam  axes  and 
receiver  beam  axis  were  constant  throughout  the  data  run. 
From  the  collection  of  scattered  sum  frequency  intensity 
spectra  (from  the  crossed  beam  experiment)  measured  across 
the  jet,  Korman  and  Rife  discovered  that  the  rms  radial 
component  of  turbulent  velocity  is  directly  proportional  to  the 
rms  acoustic  pressure  of  the  scattered  sum  frequency.  This 
relationship  demonstrated  the  apparatus's  ability  to  spatially 
resolve  components  of  the  turbulence  along  the  line  of  forward 
scattering  (or  the  line  of  bisection  of  the  primary  beam  axes). 

In  addition,  Korman  and  Rife  found  that  the  skewness  and 


kurtosis  of  the  scattered  spectra  matched  the  expected 
properties  of  skewness  and  kurtosis  of  the  turbulent  radial 
velocity  correlations  distributed  across  the  jet.  They  verified 
that  the  spectral  shape  must  be  related  to  the  probability 
density  function  for  the  turbulent  velocities.  In  contrast,  the 
scattered  spectra  for  the  "conventional"  scattering  of  a  single 
primary  beam  experiment  showed  much  less  correlation  with 
the  turbulent  flow.  It  lacked  good  spatial  resolution  and  could 
not  predict  either  radial  or  axial  turbulent  flow  profiles  across 
the  jet.  The  usefulness  of  nonlinear  crossed  beam  scattering 
over  "conventional"  single  beam  scattering  is  now  apparent. 

The  successes  of  these  earlier  works  have  motivated  this 
focused  beam  experiment. 

In  our  experiment,  the  intensity  spectra  are  measured  as 
a  function  of  angle  at  many  scan  positions  across  the  width  of 
the  turbulent  jet.  This  data  will  be  used  to  determine  the 
nonlinearly  scattered  sum  frequency  intensity  spectrum  versus 
angle  as  well  as  the  statistical  nature  of  the  turbulent  velocity 
components  (since  the  turbulent  velocity  distribution  in  the  jet 
is  anisotropic).  Measurements  of  the  Doppler  shift,  variance, 
skewness,  and  kurtosis  of  the  scattered  spectra  are  used  to 
predict  the  statistical  turbulent  velocity  information  involving 
the  axial  and  radial  direction  of  flow.  The  exact  relations 
between  the  nonlinearly  scattered  spectral  shapes  and  the 
turbulent  velocities  will  be  fully  developed  in  the  theory 
section  of  this  paper. 
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EXPERIMENTAL  SETUP 

A.  Scattering  geometry 

The  geometry  of  the  scattering  problem  involving  two 
ultrasonic  focused  beams  and  turbulence  is  chosen  to  reflect  a 
high  degree  of  symmetry.  The  mutually  perpendicular  focused 
beams  are  generated  from  two  individual  transducer  units  of 
identical  construction  Each  transducer  unit  has  a  focal  length 
of  14  cm.  The  beams  are  aligned  to  be  mutually  perpendicular 
and  to  overlap  at  a  common  focal  point  (Fig.  3). 


FIG.  3.  Arrangement  of  focused  transducers 

Both  transmitting  transducers  generate  continuous  wave 
focused  beams  of  frequencies  foi  =  2.0  MHz  and  f02  =  2.1  MHz. 
respectively.  Focused  beams  are  produced  because  the  face  of 
each  transducer  forms  a  spherical  concave  surface.  The 
diameter  of  each  surface  is  2.54  cm  .  The  transducers  are 
designed  to  resonate  at  2  MHz.  However,  the  low  Q  factor 
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allows  them  to  be  driven  off  resonance  and  still  generate  large 
pressure  amplitudes.  Each  transducer  unit  is  submerged  in  the 
tank  and  individually  attached  to  a  radius  arm  that  is 
mechanically  suspended  above  the  tank  on  a  rotary  table.  A 

third  receiving  transducer  unit  is  placed  in  the  tank  outside  the 
interaction  region.  The  receiver  is  an  unfocused  circular  plane 
array  transducer  of  diameter  =  2.54  cm.  This  unit  is  designed 
and  tuned  for  maximum  detection  response  at  a  frequency  of 
f0+  =  foi  +  fo2  =  4.1  MHz.  The  low  Q  factor  of  the  receiver  unit 
allows  the  detection  of  a  small  range  of  frequencies  around  4.1 
MHz  with  a  flat  frequency  response.  The  receiving  transducer 
is  not  focused  because  it  is  extremely  difficult  to  keep  three 
focal  points  in  alignment  for  long  periods  of  time.  Although  the 
receiver  is  not  focused,  it  is  directional  with  a  half  power  beam 
full  width  of  0.85  degrees.  This  directional  axis  is  arranged  to 
meet  perpendicularly  with  the  axis  of  the  submerged  circular 
jet  and  cross  the  jet  at  the  interaction  region  (Fig.  4).  The 
distance  from  the  face  of  the  receiving  transducer  to  the 
interaction  region  is  15  cm. 


nozzle 


Receiver 


23  cm 


interaction 

region 


FIG  4.  Arrangement  of  receiving  transducer  relative  to  the  jet 
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This  geometry  is  useful  because  the  resulting  scattering  of  the 
sum  frequency  components  have  their  rays  perpendicular  to 
the  jet’s  mean  flow.  Therefore,  this  component  will  not 
undergo  further  Doppler  shifting  as  it  propagates  outside  of  the 
interaction  region.  Thus,  one  does  not  have  to  worry  about  the 
effects  of  rescattering  of  the  sum  frequency  component. 

Using  this  apparatus,  two  experiments  are  performed. 

The  first  is  a  fixed  angle  scan  across  the  width  of  the  jet, 
similar  to  Rife's  work.  Here,  the  two  senders  and  receiver 
remain  at  fixed  angle  relative  to  each  other  while  they  are 
translated  across  the  jet.  These  "translation  scans",  as  depicted 
in  Fig.  5,  move  the  interaction  region  laterally  to  scan  the 
turbulence  perpendicular  to  the  mean  flow.  The  x-axis  in  Fig.  5 
is  along  the  radial  direction. 


FIG  5.  Geometry  of  transducers  and  jet  during  translational 


scans 


In  the  rectangular  coordinates  given  above,  the  overlap  region 
maintains  a  constant  position  of  z  =  0  d  and  y  =  33.9  d  while  x 
traverses  from  -  8  d  to  +8  d,  where  d  represents  the  diameter 
of  the  nozzle  exit  (0.25  inch  =  0.635  cm).  Two  experiments  are 
performed:  "forward  scattering"  with  the  two  senders  facing 
towards  the  receiver  (  0j  =  45°,  02  =  -45°  )  and  "back  scattering" 
with  the  two  senders  facing  away  from  the  receiver 
(0j  =  225°,  02  =  135°).  Fig.  6  illustrates  the  relative  transducer 
positions  for  "forward"  and  "back"  scattering. 
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FIG.  6.  Geometry  of  forward  and  back  scattering  experiments 

The  second  experiment  involves  rotating  the  two 
transmitting  transducers  in  a  plane  parallel  to  the  water's 
surface.  Throughout  the  experiment,  the  radius  arms  maintain 
the  beam  axes  at  right  angles  to  each  other  while  allowing 
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them  to  pivot  about  an  axis  perpendicular  to  the  surface  of  the 
water  and  through  the  point  of  overlap.  The  receiver  remains 
stationary  at  one  scan  position  along  the  x-axis  as  the 
scattering  angle  is  changed.  A  full  360®  of  scattering 
measurements  is  completed  in  3®  increments.  After  the 
apparatus  has  completed  a  360®  rotation,  (which  we  shall  call  a 
complete  data  run)  all  three  transducers  are  translated  0.2 
inches  along  die  x-axis  using  the  mechanical  translating 
carriage  located  above  the  tank.  Then  another  data  run  is 
performed.  We  completed  twelve  successful  runs  covering  0.0 
to  2.0  inches  along  the  x-axis  in  steps  of  0.2  inches.  Note  that 
gaps  exist  in  the  data  at  certain  angular  sectors.  These  gaps  are 
caused  because  the  transducers  would  be  placed  in  the  jet  if  a 
scan  was  performed  at  that  angle.  This  would  interfere  with 
the  jet  and  knock  the  transducers  out  of  alignment  so  these 
sectors  were  skipped.  Each  data  run  took  two  or  three  days  to 
complete.  This  is  due  partly  because  at  each  angle  it  took  45 
seconds  to  electronically  sweep  and  detect  one  spectrum.  It 
was  necessary  to  sweep  20  times  to  get  one  average  spectrum. 

B.  Transmitting  and  Receiving  Electronics 

Measurements  of  the  Doppler  shift  in  the  sum  frequency 
spectrum  require  an  absolute  calibration  of  frequency  of  10  Hz 
in  4.1  MHz.  To  insure  that  no  drift  occurs  in  our  primary 
signals  (that  will  raise  questions  upon  the  accuracy  of  the 
Doppler  measurements),  it  is  necessary  to  require  excellent 
frequency  stability.  Two  entirely  separate  but  nearly  identical 
crystal  oscillator  circuits  are  designed  to  generate  the  two 
primary  electronic  frequencies  foi  (  =  2.006944  MHz)  and  fo2 
(  =  2.100609  MHz)  with  high  stability.  Fig.  7  presents  an 
overview  of  the  transmitting  electronics  in  block  format. 


FIG  7.  Transmitting  electronics 

A  buffer  amplifier  protects  the  crystal  oscillator  from  voltage 
fluctuations  which  could  affect  the  amplitude  and  frequency  of 
the  oscillator's  signal.  Next,  the  signal  is  fed  through  50  ohm 
five  stage  Butterworth  low  and  high  pass  filters  (2.25MHz  and 
1.75  MHz  cut-off  frequencies,  respectively)  to  insure  the 
signals  are  pure  sinusoids  with  virtually  no  harmonic 
distortion.  Finally,  a  50  ohm,  broad  band,  100W  radio 
frequency  power  amplifier  increases  the  signal's  amplitude. 

An  L-C  matching  and  tuning  circuit  is  used  to  couple  this  power 
into  each  transmitting  transducer  unit.  Thus,  the  sound 
generated  is  a  pure  tone  at  a  frequency  equal  to  fo i  or  fo2-  The 
electronic  signal  amplitude  is  measured  to  be  60  volts  peak  to 
peak  across  the  transducer  input  cable. 

The  receiving  transducer  and  electronic  instrumentation 
detect  the  relatively  weak  sum  frequency  pressure  while 
filtering  out  the  intense  primary  frequencies  fo i  and  fo2  that 
can  also  insonify  the  receiving  transducer.  A  block  diagram  of 
the  receiving  electronics  is  presented  in  Fig.  8. 


FIG.  8.  Receiving  Electronics 

A  50  ohm  BNC  cable  connects  the  4.1  MHz  receiver  to  a  high 
pass  filter.  This  5  stage  50  ohm  Chebyshev  passive  high  pass 
filter  has  a  3  dB  cutoff  frequency  at  3.75  MHz.  The  filter 
insures  that  no  2.0  or  2.1  MHz  signals  are  present  in  the 
amplifier  or  spectrum  analyzer  electronics.  If  these  low 
frequency  signals  are  allowed  to  exist,  the  spectrum  analyzer 
circuitry  will  electronically  mix  them  nonlinearly,  which  will 
produce  a  sum  frequency  on  its  display.  This  electronically 
generated  sum  frequency  is  highly  undesirable  since  it  would 
deteriorate  the  ability  to  measure  the  acoustically  generated 
sum  frequency  component. 

Since  this  laboratory  facility  was  not  shielded  from 
outside  radio  frequency  transmissions,  all  electronic  devices 
are  placed  in  metal  shielded  boxes,  and  all  connections  between 
components  are  made  with  50  ohm  coaxial  cables.  This 
precaution  insures  that  the  signals  received  by  the  spectrum 
analyzer  are  generated  acoustically  at  the  receiving  transducer. 
After  the  high  pass  filter,  a  low  noise,  linear,  radio  frequency 
amplifier  amplifies  the  sound  signal  and  outputs  it  to  a 
Tektronix  7L5  Sweep  Spectrum  Analyzer.  Here,  the  frequency 
components  of  the  incoming  signal  are  decomposed  using  a  30 
Hz  resolution  sweep  filter.  The  voltage  amplitude  at  each 
frequency  is  digitized  by  an  IQ400  Digital  Oscilloscope  which 


then  averages  the  results  of  twenty  signal  sweeps  and  stores 
the  resulting  spectra  to  disk  on  a  Macintosh  Plus  computer  for 
later  analysis.  The  IQ400  digitized  at  12  bits  over  a  +  200  mV 
range  for  1024  points,  with  a  sample  time  interval  of  50  msec.§ 

C.  Turbulence 

A  submerged  circular  water  jet  generates  the  turbulence 
in  our  experiment.  This  jet  is  formed  from  a  1/4  inch  (0.635 
cm)  diameter  nozzle  and  powered  by  a  1/8  hp  centrifugal 
water  pump.  The  pump  creates  4.15  psi  of  pressure  which  was 
used  to  predict  a  nozzle  exit  velocity  of  Uo  =  7.1  m/s.  Care  is 
taken  to  insure  that  there  is  a  high  degree  of  laminar  flow  at 
the  nozzle  exit.  Therefore,  turbulent  shear  flow  is  generated  by 
fluid  flow  mechanisms  that  are  dominated  by  entrainment, 
turbulent  mixing,  vortex  stretching,  and  nonlinear  diffusion 
mechanisms.  Care  was  taken  so  that  turbulent  flow  does  not 
develop  from  turbulence  that  already  exists  at  the  jet  exit. 

Such  "pre-jet"  turbulence  can  be  created  from  turbulent  flows 
in  the  pump  hoses  and  fittings  leading  to  the  nozzle.  A  "clean 
jet"  or  laminar  jet  at  the  nozzle  exit  is  accomplished  by 
introducing  a  plenum  chamber  prior  to  the  exit.  Fig.  9  shows 
the  location  of  the  plenum  and  other  parts  of  the  turbulence 
generation  apparatus. 


5  In  designing  the  experiment,  care  was  taken  to  choose  a 
sum  frequency  component  that  would  not  be  close  to  the  high 
powered  radio  frequencies  generated  at  the  North  Severn 
Naval  Station,  Annapolis,  MD.  This  fact  had  forced  a  design 
change  from  a  planned  4.0  MHz  sum  frequency  to  4.1  MHz. 


FIG.  9.  Turbulence  generation  apparatus 


The  plenum  mitigates  any  turbulent  motion  by  allowing 
the  flow  cross-section  to  expand  over  an  extended  region. 

From  the  equation  of  continuity  for  incompressible  fluid  flow, 
the  ratios  of  the  cross-sectional  areas  at  two  points  in  a  flow 
must  be  inversely  proportional  to  the  velocities  at  those  same 
two  points.  The  cross-sectional  area  of  the  plenum  is  36  sq. 
inches.  The  one  inch  diameter  hose  leading  into  it  has  a  cross- 
section  of  0.78  sq.  inches.  Therefore,  the  velocity  of  the  water 
must  decrease  by  a  factor  of  about  46.  The  turbulent  motion  in 
the  plenum  (a  small  fluctuation  of  the  mean  flow)  will  reduce 
by  the  same  factor. 

The  plenum  connects  to  a  conical  section  (Fig.  10)  that 
reduces  the  flow  from  a  6  inch  diameter  to  a  3  inch  diameter 
over  a  length  of  6.25  inches. 
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conical  section  machined 

noz2le 


FIG.  10.  Nozzle  construction 

This  conical  section  connects  to  a  high  strength  epoxy  nozzle 
(Fig.  10)  that  gradually  reduces  the  3  inch  diameter  flow  down 
to  a  1/4  inch  over  a  1.5  inch  distance.  This  nozzle  is  poured 
and  set  into  a  machined  aluminum  mold  that  is  manufactured 
and  shaped  on  a  computer  driven  lathe  cutter.  The  machining 
was  performed  in  the  Rickover  Hall  Machine  Shop  under  the 
supervision  of  Mr.  Carl  Owen.  The  nozzle  shape  is  designed  to 
keep  the  flow  highly  laminar  in  the  nozzle  throat.  The 
turbulence  from  this  jet  (laminar  at  the  exit)  can  now  be 
compared  with  other  jets  that  are  documented  in  the  literature. 

D.  Mechanical  and  Computer  Control 

The  experiment  is  performed  with  the  assistance  of 
computer  control  and  mechanical  automation  because  of  the 
length  of  time  required  to  complete  a  data  run.  Translational 
runs  (at  a  fixed  angle)  take  twenty  four  hours  each,  while 
rotational  scans  (at  a  fixed  scan  position)  require  several  days. 
An  Apple  He  is  chosen  for  the  computer  control  because  of  its 
ease  in  programming  and  for  its  ability  to  interface  with 
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several  very  different  machines.  The  Apple  lie  controls  a 
stepper  motor  to  position  the  transducers  at  the  appropriate 
scattering  angle  for  each  data  sweep,  while  for  the  translational 
runs.  Another  stepper  motor  is  controlled  by  the  Apple  lie  to 
move  a  piston  that  in  turn  pushes  a  carriage  with  grooved 
wheels  mounted  on  two  inverted  "V"  tracks  that  are  supported 
above  the  tank.  The  orientation  of  the  carriage  and  tracks  with 
respect  to  the  jet  is  shown  in  Fig.  11. 


Rails 


FIG.  11.  Translational  carriage  orientation 

Since  the  tracks  are  parallel  to  the  x-axis,  the  carriage  moves 
all  three  transducers  along  the  width  of  the  jet.  In  the 
rotational  mode,  a  stepper  motor  and  a  reduction  gear  box  turn 
the  crankshaft  on  a  rotary  table.  The  radius  arms  of  the 
sending  transducers  are  mounted  on  this  table.  The  two 
senders  rotate  in  a  plane  parallel  to  the  water's  surface. 

Communication  with  the  spectrum  analyzer  involves 
triggering  the  spectrum  analyzer’s  single  sweep  control.  This  is 


accomplished  by  dedicating  a  digital  switch  to  trigger  the  7L5 
at  the  appropriate  time  determined  by  an  Apple  program.  A 
ramp  signal  from  the  7L5  sweep  time  base  then  triggers  the 
IQ400  to  trigger  its  acquisition  mode.  The  Apple  lie  also 
controls  the  IQ400  by  an  IEEE-488  GPIB  interface  bus.  This 
link  allows  the  Apple  to  control  the  IQ400  to  signal  average  the 
incoming  spectra,  direct  the  IQ400  to  save  each  spectrum  to  a 
Macintosh  Plus  computer  for  future  analysis,  and  keep  track  of 
scan  position,  rotation  angle  and  the  labeling  of  the  data.  The 
computer  programs  and  associated  documentation  are 
presented  in  appendix  C. 


THEORETICAL  DISCUSSION 
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A.  Overview 

The  theory  involving  the  nonlinear  interaction  of  crossed 
ultrasonic  beams  is  greatly  simplified  if  the  interaction 
involves  incoming  plane  waves.  In  order  to  use  Korman  and 
Beyer's  theory5,  we  must  model  the  focused  incoming  beams  as 
plane  waves  inside  the  tightly  defined  interaction  region  where 
the  two  beams  overlap.  The  focused  transducers  do  not 
actually  emit  plane  waves;  however,  since  the  curvature  of  the 
wave  fronts  is  small  compared  to  the  size  of  the  overlap  region, 
the  plane-wave  approximation  is  acceptable.  Thus,  we  are 
justified  in  using  Korman  and  Beyer's  plane  wave  theoryS  for 
our  experiment. 

Of  particular  importance  to  the  experiment  is  the  Doppler 
shift  equation,  tod  =  •  V,  which  relates  the  magnitude  of  the 

Doppler  shift  to  the  sum  frequency  component's  wavenumber 
and  the  velocity  of  the  turbulence.  We  will  now  discuss  the 
Doppler  phenomenon  in  general  so  that  the  reader  can  become 
familiar  with  its  use  in  the  crossed  beam  experiment. 

A  Doppler  shift  is  a  change  of  frequency  due  to  the 
relative  motion  between  a  sound  source  and  a  receiver.  A 
common  example  is  the  rising  pitch  of  a  train  whistle  as  it 
approaches  a  stationary  observer  at  a  train  station.  In  the  case 
where  the  train  is  approaching,  the  train’s  whistle  closes  in  on 
the  sound  waves  in  front  of  it.  The  wavelengths  of  sound  then 
become  shorter.  This  increases  the  sound’s  frequency  which 
you  hear  as  a  higher  pitch.  In  our  case,  the  sources  are  tiny 
packets  of  sound  wavelets  formed  in  the  overlap  region  by  the 
interaction  of  the  crossed  sound  beams  and  the  turbulent 
scattering  sites.  These  packets  of  sound  sources  (or  phonons) 
are  jostled  around  by  the  action  of  the  turbulence  which 
creates  a  relative  motion  between  the  wavelets  and  the 
stationary  4.1  MHz  receiver.  Thus,  although  the  "phonons"  are 


theoretically  emitting  a  pure  4.1  MHz  tone,  their  relative 
motion  toward  or  away  from  the  receiver  causes  the  receiver 
to  detect  a  Doppler  shifted  frequency  that  can  be  higher  or 
lower  than  the  4.1  MHz  tone.  The  Doppler  shifted  frequency  is 

IC  •  V 

given  by  f0+  +  f<j  where  fd  =~  =  ~ - . 

ZK  L  K 

This  equation  for  Doppler  shifts  is  used  to  relate  the  shape  of  a 
scattered  sum  frequency  spectra  to  the  time  averaged 
turbulent  velocity  quantities. 

Unlike  the  train  whistle  which  has  only  one  shift  in 
frequency,  the  scattered  sum  frequency  has  several  shifts. 
Correctly  speaking,  it  has  a  band  of  shifts.  Each  scattered 
wavelet  is  caused  by  a  turbulent  eddy  in  the  overlap  region 
that  has  a  characteristic  velocity.  Since  the  turbulent  eddies  or 
scattering  sites  that  generate  different  waves  varies  spatially 
over  the  interaction  region  and  with  a  velocity  that  can  change 
with  time,  each  wavelet  is  jostled  by  a  different  amount 
depending  upon  the  local  motion  of  the  turbulent  jet.  This 
causes  several  different  tones  to  be  detected  at  the  receiver. 

All  tones  are  around  fo+  but  some  are  above  it  (the  turbulent 
eddy  moves  toward  the  receiver)  and  others  are  below  (the 
turbulent  eddy  moves  away  from  the  receiver).  The  net  result 
is  to  smear  the  pure  tone  fo+  over  a  band  of  frequencies 
centered  around  fo+.  Therefore,  the  net  result  is  that  the 
scattered  sum  frequency  component  is  detected  as  a  scattered 
spectrum  of  frequencies.  This  spectrum  is  called  the 
nonlinearly  scattered  intensity  spectrum. 

The  shape  of  the  scattered  spectra  is  quantified  by 
calculating  the  average  frequency  <  f  >,  variance  o2  ,  skewness 
S,  and  kurtosis  K  of  the  scattered  spectra.  In  our  research 
group  these  are  called  the  shape  factors.  These  spectral 
quantities,  (<  f  >,  o2,  S,  and  K)  each  describe  a  specific 
contribution  to  the  broadened  scattered  sum  frequency 
spectrum  which  can  be  related  to  the  vector  turbulent 
velocities  that  are  the  cause  of  the  Doppler  shifts.  The  average 
(or  most  probable)  frequency  in  the  scattered  spectrum  is  due 


to  the  mean  or  average  velocity  of  the  turbulent  eddies  in  the 
overlap  region.  This  frequency  is  given  by  <  f  >  =  <fo+  +  /  ~^> 

=  fo++<lO  V>/2n,  where  V  is  the  instantaneous  velocity  of  the 

scattering  eddies.  If  V=V0  +  v,  where  V0  is  the  average 

velocity  and  V  is  the  fluctuating  portion,  then  the  time  average 
Doppler  shift  <  f  >  can  be  expressed  by  <  f  >  -  fo+  which  is  <fd>  =< 
f  >  .  f^  =  (K>  <V>)/2ji. 


FIG.  12.  Illustration  of  Doppler  shift  and  variance 

In  Fig.  12,  the  broadened  curve  represents  the  measured 
scattered  intensity  spectrum  versus  frequency.  The  spike  at 
fo+  is  the  pure  tone  of  the  sum  frequency  which  has  been 
smeared  out  to  give  the  observed  spectrum.  Define  <  f  >  to  be 
the  average  frequency  of  the  spectrum.  It  is  determined  from 
the  mathematical  expression 

J  Kf)  fdf 


<{>  = 


J  KO  df 


One  can  compute  the  variance  of  the  spectrum.  The 
standard  deviation,  o,  is  the  rms  frequency  and  is  a  measure  of 
the  width  of  the  spectrum.  Mathematically,  o:  is  given  by 


o2  = 


J 


1(f)  (f  -  <  f  >)2  df 


J 


1(f)  df 


The  magnitude  of  the  turbulent  velocity  f  m  an 

interaction  region  determines  the  broadening  jf  the  sum 
frequency  spectra.  If  the  velocity  fi  equations  in  the 
interaction  region  are  large,  th^n  the  distribution  of  Doppler 
shifts  "heard"  and  the  spectral  broadening  will  be  large  as  in 
the  diagram  at  the  left  in  Fig.  12.  Conversely,  if  there  are  only 
small  velocity  fluctuations,  then  there  will  be  only  a  small 
distribution  of  Doppler  shifts  and  the  spectra  will  be  narrow 
which  is  shown  by  the  diagram  at  the  right  in  Fig.  12. 

Skewness  indicates  a  bias  in  the  spectra  towards  one  side 
of  the  frequency  spectrum  or  the  other.  Intensity  spectra 
exhibiting  positive  and  negative  skewness  are  shown  in  Fig.  13. 


FIG.  13.  Comparison  of  positive  and  negative  skewness 

Skewness,  S,  is  determined  mathematically  from  the  third 
moment  of  the  spectrum.  It  is  given  by 


2  5 

1(f)  (f  -  <  f  >)3  df 


For  a  gaussian  or  any  symmetric  distribution  function,  the 
skewness  is  always  equal  to  zero,  indicating  no  bias.  If  the 
turbulent  velocity  fluctuations  spend  more  time  moving  toward 
the  receiver  than  away  from  it,  then  the  intensity  spectrum 
will  exhibit  a  positive  skewness,  S.  In  contrast,  a  negative 
skewness  means  that  the  velocity  fluctuations  spend  more  time 
moving  away  from  the  receiver  than  toward  it. 

The  last  spectral  shape  property  is  called  the  Kurtosis 
which  describes  the  intermittency  or  bursting  phenomenon  of  a 
probability  distribution. 


FIG.  14.  Comparison  of  distributions  exhibiting  large  and  small 

kurtosis 

Define  the  kurtosis,  K.  to  be  given  by 
f  1(f)  (f  -  <  f  >)4  df 


In  Fig.  14,  the  graph  at  left  exhibits  a  small  kurtosis.  This 
distribution  represents  a  turbulent  fluctuation  that  is  steady 
and  does  not  change  much  in  time.  In  contrast,  the  graph  on 
the  right  in  Fig.  14  is  related  to  a  distribution  with  a  very  large 
kurtosis.  This  would  describe  a  particle  that  has  at  first  a 


period  of  steady  fluctuations,  then  suddenly  exhibits  large 
bursts  of  velocity  changes  which  later  go  back  to  being  steady 
fluctuations.  For  a  Gaussian  distribution,  the  kurtosis  is  equal 
to  three,  which  is  often  used  as  a  comparison  to  other 
distributions. 

With  these  descriptions  of  the  spectral  distributions,  we 
can  discuss  the  nonlinear  crossed  beam  spectral  theory.  From 
our  initial  hypothesis  that  the  focused  beams  will  create  locally 
plane  waves  within  the  overlap  region,  one  can  make  use  of 
Korman  and  Beyer's  theory5  with  plane  waves.  From  their 
nonlinear  scattering  theory,  they  derived  an  equation  to  define 
the  Doppler  shifts  of  the  scattered  sum  frequency  component, 
tod=  K+  •  V.  Using  this  equation,  one  can  extract  values  for  the 
turbulent  velocity  fluctuations  (as  statistical  quantities)  from 
the  experimentally  determined  scattered  sum  frequency 
spectrum  measured  as  a  function  scattering  angle.  To  describe 
this  method,  one  first  defines  these  symbols: 

cod  =  Doppler  shift  of  the  sum  frequency  in  rad/sec 
K+  =  wavenumber  vector  of  the  scattered  sum 
frequency  component 

=  k[  +  k2 

Kj  =  scattered  wavenumber  vector  corresponding  to 
scattering  from  a  conventionally  scattered  single 
sound  beam  originating  from  Transducer  #1 
__  =  koi(n  -  "Ol) 

K2  =  similar  definition  for  the  sound  beam  emitted 
by  Transducer  #2 

koi,  k<)2  =  magnitude  of  the  wavenumber  from 
transducer  #1,  transducer  #2 
=  ci>oi/c,  CO02/C 

cdoi,coo2  =  an8u*ar  frequency  (rad/sec)  emitted  from 
transducer  #1,  transducer  #2 

c  =  speed  of  sound  in  water  (1482  m/s  at  20°C) 
n  =  unit  vector  pointing  towards  the  receiver 
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noi  =  unit  vector  pointing  along  original  path  of  beam#l 
n02  =  unit  vector  pointing  along  original  path  of  beam#2. 

Fig.  15  illustrates  the  relationships  between  n,  noi,  and  no2- 


FIG.  15.  Incident  and  scattered  sound  beam  geometry 

Both  transmitting  transducers  are  free  to  rotate  in  the 
plane  described  in  Fig.  16.  The  angles,  0j  and  02,  are  used  to 
define  their  angular  positions  with  respect  to  the  receiver  axis 
and  each  of  the  respective  transmitting  axes.  Since  the  relative 
angle  between  the  two  transmitting  transducers  is  always 
maintained  at  90\  a  single  scattering  angle,  0,,  can  be  used  to 
describe  the  transducer  positions.  Express  the  scattering  angles 
0!  and  02  in  terms  of  the  symmetry  angle  0,  by 
0,  =  0.  +  45* , 

02  =  0.  -  45*  . 

This  symmetry  is  shown  in  Fig.  16. 
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FIG.  16.  Illustration  of  the  angular  relationships 

It  is  important  to  recognize  that  two  rectangular 
coordinate  systems  must  be  developed  to  completely  define  the 
scattering  geometry  (Fig.  17).  The  first  matches  the  usual  polar 
coordinate  system  set  up  by  the  orientation  of  the  transducers. 
The  second  defines  the  normal  coordinate  system  from  the  jet's 
point  of  view. 
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FIG.  17.  Comparison  of  normal  coordinate  systems 

Notice  the  inversion  of  the  x-y  plane  between  the  two  systems. 
Measurements  in  the  x-y  plane  of  one  system  will  have  the 
opposite  sign  in  the  other.  This  forces  a  negative  sign  to  the 
equation  for  Doppler  shifts  to  maintain  consistency  between 
the  systems,  so  tod  =  -(K+  •  v) . 
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To  model  the  turbulent  flow,  let  V  =  (Vqx  +  Vx)i  +  (Voy 

+Vy)j 

where 

V0x  =  mean  velocity  in  the  x  direction 

V0y  =  mean  velocity  in  the  y  direction 

Vx  =  turbulent  velocity  fluctuations  in  the  x  direction 

Vy  =  turbulent  velocity  fluctuations  in  the  y  direction. 

The  turbulent  velocity  fluctuations  average  to  zero.  This  is 
expressed  by  <VX>  =  0,  and  <Vy>  =  0,  where  the  brackets  around 
a  variable  denote  an  average  of  that  variable  over  time. 

Now  the  definitions  of  <  f  >,  o2  ,  S,  and  K  can  be  related  to 
four  different  spectral  moments  which  are  constructed  from  cod- 
It  will  be  shown  later  that  these  four  moments  will  yield  linear 
combinations  of  turbulent  velocity  moments.  These  moments 
are  often  referred  to  as  correlation  coefficients  of  the 
turbulence. 

The  theory  is  now  outlined  below: 


=  -(*♦  «_V) 

=  -  (Kj  +  K2)  •  V 

*  koi(n  -  not) 

*  ko2(n  -  n02) 

=  (Vox  +  Vx)i  +  (Voy  +  Vy)j 

not  =  i  cos(6j)  +j  sin(0!) 
n  o  2  =  i  cos(02)  +  j  sin(02) 


from  Korman  and  Beyer5  (1) 

definition 

definition 

turbulent  eddy's  velocity 
vector 

unit  vectors  in  the  direction 
of  each  incident,  primary 
wave 

(2) 

(3) 

(4) 


=  koi[(l-cos(0,))  i-  sin(0!>j] 

^2  kO2[(l-cos(02))  i-  sinCO^j] 

Kj  +  K2  =  [koi(l-cos(0j)  +  ko2(l-cos(02)]  i 
-  [koisin(0,)  +  kozsinC©^]  j  . 


CDd 

K2 

V 


We  can  write  this  expression  in  terms  of  0*  by  using  this 
transformation: 


e2=e.-45~ 

fC  =  ko+(l  -  e.  -Tsin  0.)  i  +  ko+(  -^«in  6.  -  Y^cos  0.)  j  (5a) 
where 

ko+  =  koi  +  ko2 

ko-  =  koi  -  ko2 

y  =  k0.  /  k0+ 

Let 

Kx  =  k0+(l  -  ^cos  0.  -  ysin  0.)  (5b] 

Ky  =  ko+(  -  ^sin  0.  -  V^cos  0.)  ,  then  (5c) 

cod  =  -[Kxi  +  Ky j]  •  [(Vox  +  Vx)i  +  (Voy  +  Vy)j]. 

By  substitution, 

cod  =  -[KX(V0X  +  Vx)  +  Ky(Voy  +  Vy)]  .  (6) 

Since  the  instantaneous  values  of  the  Doppler  shift 
depend  upon  the  instantaneous  velocities  of  the  turbulence, 
which  are  constantly  fluctuating,  we  are  more  interested  in 
time  averaged  values, 

<  <*>  =  -[Kx(Vox  +  <VX>)  +  Ky(Voy  +  <Vy>)] . 

Since  Kx,  Vox,  Ky,  and  Voy  are  all  constant  with  respect  to  time, 

<a*>  =  -(KxVox+KyVoy).  (7) 

Now  we  relate  the  average  Doppler  shift  to  the  average 
frequency. 


<f>  =f0+-<fd> 


(8) 


<f>  =  fo+  -  <  G>d  >/2jc  ; 
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which  yields  a  relation  from  which  we  can  extract  the  mean 

velocity  of  the  jet  using  our  scattered  spectra 

<f>-fot  =  (K,V„  +KyVoy)/2*  .  (9) 

To  set  up  the  second,  third,  and  fourth  moments,  define 
the  quantity  o>d  -<  oid  >  below: 

©d-<tOd>  =-(KxVx+KyVy).  (10) 

The  second  moment  denoted  by  a2  =  <(cad  -<  cod  >)2>. 
becomes 

<(©d  -<  ©d  >)2>  =  <(KXVX  +  Ky  Vy)2> 

=  Kx2< Vx2>  +  2KxKy< Vx  Vy>  +  Ky2<Vy h> .  (11) 

Similarly,  the  third  moment  =  <(©d  -<  ©d  >)3>  becomes 
<(©d  -<  oid  >)3>  =  <(KXVX  +  KyVy)3> 

=  -  KX3<VX3>  -  3Kx2Ky< VX2V y>  -  3KxKy2<VxVy2> 

-  Ky3<Vy3> .  (12) 

Using  this  relation,  the  skewness  is  defined  to  be 

S  =  <(©d  -<  a>d  >)3>  /  a3. 

Finally,  the  fourth  moment  =  <(©d  -<  “d  >)4>  becomes 

<(©d  -<  oid  >)*>  =  <(KXVX  +  KyVy)4> 

=  KX4<VX4>  +  4KX3  Ky<  Vx3  V  y>  +  6Kx2Ky2<Vx2Vy2>  + 

4KxKy3<VxVy3>  +  Ky4<Vy4>.  (13) 

Kurtosis  is  defined  as: 


K  =  <(oid  -<  oid  >)4>  /  a4  . 
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From  these  relations  of  the  average  Doppler  shifts  (which 
come  from  statistical  properties  of  the  broadened  sum 
frequency  spectra),  one  can  find  Vox,Voy,  <VX2>,  <VxVy>,  <Vy2>, 
<VX3>  ,  <Vx2Vy>,  <VXVy2>,  CVy^, 

<VX4>,  <Vx3Vy>,  <Vx2Vy2>,  <VxVy3>,  and  <Vy4> 
of  the  turbulent  flow. 

B.  Isotropic  turbulence 

If  one  assumes  that  the  turbulence  is  isotropic  with 
respect  to  the  randomly  fluctuating  velocity  component,  then 
the  magnitudes  of  the  fluctuations  in  the  radial  (x)  and  axial  (y) 
directions  are  equal.  This  isotropic  model  allows  the  equations 
for  a2,  S,  and  K  to  be  simplified.  Although  this  is  an 
oversimplification,  the  resulting  relationships  between  spectral 
moments  and  turbulent  velocities  are  easy  to  visualize. 

Let  V  equal  the  magnitude  of  the  velocity  fluctuations 
with  components  Vx  and  Vy  given  by 
Vx  =  Vcos  a 

Vy  =  Vsin  a  , 

where  V  and  a  are  functions  of  time. 


Here,  one  imposes  the  condition  of  isotropy  by  stating  that  V 
and  a  are  statistically  independent  random  variables.  Using 
this  model,  the  variance,  a2,  becomes 


a2  =  Kx2<V2cos  a2>  +  2KxKy<V2cos  asin  a>  +  Ky^V^in  a2> . 


Since  <cos  a2>  =1/2,  <cos  asin  a>  =  0,  and  <sin  a2>  =  1/2, 


2^  - 


=  <V2>[Kx2  (1/2)  +  0  +  Ky2(l/2)] 
=  (1/2)<V2>[Kx2+  Ky2]  . 


(14) 


Since  Kx2  +  Ky2  is  constant  at  a  fixed  scattering  angle,  a2  i  s 
proportional  to  the  square  of  the  turbulent  velocity. 

The  expression  for  skewness  becomes 


S  =  (-  Kj^cV^  <cos  cx'b*  -  3Kx2Ky<V:J>  <cos  a2sin  a  >  - 
SKxKy^V^  <cos  asin  a  2>  -  Ky3<V3>  <sin  a3>)  /  a3  . 


Since  <cos  ah>,  <cos  a2  sin  a>,  <cos  asin  ab>,  <sin  a3>  all  =  0  , 


=  0 


The  expression  for  kurtosis  becomes 


K  =  Kx4<V4>  <cos  0'S  +  4Kx3Kv<V4>  <cos  a3  sin  a>+ 


bK^Ky^V^  <cos  a2  sin  a2>  +  4KxKy3<V4>  <cos  asin  a^  + 
<sin  a4>  /  c4. 


Since  <cos  a4>  =  3/8,  <cos  a3sin  a>  =  0,  <cos  a2sin  a2>  =  1/8, 
<cos  asin  a^  =  0,  <sin  a4>  =  3/8, 


(15) 


K 


=  <V4>[Kx4(3/8)  +  0  +  6Kx2Ky2(l/8)  +  0  +  Ky4(3/8)]  /  o4 
=  <V4>  (3/2)  [( 1  /4)(KX4  +  2Kx2Kv2  +  Kv4)] 

<V2>2  [(1/4)(KX4  +  2Kx2Ky2  +  Ky4)] 

=  (3/2)  Kturb . 


(16) 


Here,  Kturb  is  the  kurtosis  of  the  isotropic  turbulent  velocity 
fluctuations  at  the  current  scan  point. 


C.  Anisotropic  turbulence 


The  actual  turbulent  flow  is  anisotropic.  Therefore,  to 
find  the  turbulent  velocities,  one  must  use  curve  fitting 
techniques  and  apply  the  theoretical  model  of  nonlinear 
scattering.  Two  methods  are  used  to  collect  nonlinear  scattered 
data.  Therefore,  there  are  two  separate  procedures  used  to 
develop  curve  fits.  Once  the  theoretical  model  is  fitted  to  the 
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actual  shape  factors  of  the  scattered  sum  frequency  spectra, 
then  the  turbulent  velocity  correlations  can  be  predicted.  Each 
method  is  described  in  the  sections  that  follow. 


APPLICATION  OF  THE  SPECTRAL  THEORY  TO  THE 
TRANSLATIONAL  CASE  OF  CROSSED  BEAM 
SCATTERING  (e.  FIXED) 


The  first  experiment  applies  the  spectral  theory  to  the 
case  of  translating  the  crossed  beams  across  the  jet  at  a  fixed 
angle.  Two  types  of  translational  runs  are  attempted,  "forward 
scattering"  in  which  0,  =  0*  and  "back  scattering"  where  0.  = 
180*.  By  solving  the  following  system  of  equations  from  the 
Doppler  shift  relation,  a  closed  form  solution  is  derived  to  find 
the  mean  axial  and  radial  velocities. 

Modifying  Eq.  7  from  the  spectral  theory  yields 

«*>  =  -[Kx(0.)Vox  +Ky(0,)Voy]  , 

which  makes  explicit  the  dependence  of  Kx  and  Ky  on  0.. 

Substitute  0.  =  0*  and  0.  =  180*  into  Eq.  7  along  with  the 
experimentally  determined  average  Doppler  shifts  to  create  a 
system  of  two  equations: 

<<ad(e.  =0*)>  =  -(Kk(0*)Vok  +Ky(0')Voy)  (17) 

<<0d(e.  =  180*)>  =  -(K*(180-)Vo:l  +  Ky(180-)Voy).  (18) 


To 

determine  Vox  and  Voy,  first  find  Kx  and  Ky  from 

Eqs(5b,  5c); 

Kx(0*) 

=  ko+  -  ?  ko+ 

Ky(0*) 

* 

i 

it 

Kx(180‘) 

=  ko+  +  ^  ^0  + 

Ky(180’) 

ii 

Substitution  of  these  values  and  dividing  through  by  coo+ 
(to  make  all  terms  dimensionless)  yields 
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<  (od(6.  =  0*)>  /  co0+  =  (Vox[l  -  _  Voy[f  ]7)  /  c  and 

<  tod(0.  =  180*)>  /  o)0+  =  (vo*U  +  +  Voy[^F)  /  c  . 

Let  gf  =  <  cod  >1(0.  =0*)/  o)0+ 

gb  =  <  cod  >1(0.  =  180*)/  0)0+ . 

The  system  can  now  be  solved  for  the  mean  velocities  which 
are  given  as 

vox  =  (gf+ 8b)c  /  2  (19) 

vo,  =  ^  [gf  +  gb  (i^a)]c .  (20) 

It  is  pointed  out  that  the  experimental  measurements  of  the 
mean  axial  velocity,  Voy,  may  be  subject  to  large  errors  if  the 
beams  are  not  aligned  perfectly.  If  we  let  0ER  be  the  size  of  the 
angular  misalignment  then,  we  can  find  a  relation  that  will 
show  how  large  our  errors  are  due  to  the  misalignment. 

For  0.  =  0* 

sin(0*  +  0£r)  *  ®er 
cos(0.  +  0ER)  =  1  . 

For  0.  =  180* 

sin(0.  +  0ER)  =  -  0£R 
cos(0.  +  0ER)  =  1  . 

By  substituting  these  expressions  into  the  definitions  of  Kx(0.) 
and  Ky(0.) 

M0-)  -  ko+0  -  %  -  'i  eER) 

-  kof(l  -  f) 

*  k0.(- 


Ky(O') 
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Kx(180*)  =  k0+(l  +  ^  -  7  eER) 

=  k0+(l  +  ^ ) 

Ky(180*)  =  ko+(feER  +y&)  . 

Using  these  relations  for  Kx(0.)  and  Ky(0.)  and  solving  Eqs.  17 
and  18  similarly  to  before,  we  get  the  following  relations: 

V0x  =  (gf+gb)c/2  (19’) 

voy  =  (i/n  £[gf  +  gb  (i^)]c ,  (20') 

where  T  =  0ER  +7  . 

Note  that  small  errors  in  alignment  introduce  no  error  in  the 
measurements  of  the  mean  radial  velocity.  However, 
alignment  errors  are  a  considerable  factor  in  the  measurements 
of  the  mean  axial  velocity.  Since  the  ratio  of  ko-/ko+  =  7  = 

0.0228  for  our  crossed  beam  experiment,  an  angular  error  of 
0.0228  radians  (1.3  degrees)  in  the  alignment  of  the 
transducers  will  introduce  a  factor  of  two  error  in  the 
measured  value  of  mean  axial  velocity. 

Further  analysis  of  the  scattering  geometry  leads  to  the 
discovery  that  the  effective  beamwidth  of  the  transmitting 
transducers  will  also  contribute  to  the  error  factor  T. 


FIG.  18.  Error  due  to  effective  beamwidth 
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Since  the  scattering  occurs  over  a  range  of  angles  (Fig.  18),  an 
error  develops  because  multiple  Doppler  shifts  are  created. 

This  error  is  related  to  the  effective  distribution  of  angles  of 
the  focused  beam.  For  the  focused  transmitting  transducers, 
the  effective  distribution  of  angles  is  on  the  order  of  the  radius 
of  the  transducer  face  divided  by  the  focal  length.  For  our 
apparatus  this  ratio  is  0.5/5. 5  which  is  equal  to  0.09091 
radians  and  is  much  larger  than  7.  Therefore  an  experimentally 
determined  mean  axial  velocity  is  not  possible  from  the 
translational  data  taken  at  0*  and  180*. 

To  optimize  the  analysis,  one  should  choose  a  set  of  two 
angles  (a.,  J3.)  such  that  in  one  scan  experiment  Kx  is 
maximized,  and  in  the  other  experiment  Ky  is  maximized.  To 
find  the  extremum  of  Kx,  find  its  derivative  with  Tespect  to  0. 
and  set  it  equal  to  zero.  Therefore 

dKx  _  ~ 

dQ  -  y-ko+sin(0.  )  -  ko-cos(0.)  =  0  ,  and 
tan(0*)  =  . 

Thus,  0*  must  be  1.8*  or  181.8*.  Use  of  the  second  derivative 
test  finds  that  181.8*  is  the  maximum  of  the  function  Kx. 

To  maximize  Ky,  follow  a  similar  procedure  to  get  the 
following  relations 

— -  =  -  ^-ko+cos(0* )  +  ^ko-sin^.)  =  0  ,  and 
d0.  L  1 

tan(  0*)  =  1/  If . 

Once  again,  after  using  the  second  derivative  test,  the 
maximum  of  Ky  is  found  to  be  268.7*. 

Thus  scanning  at  the  angles  0.  =  181.8*  and  268.7*  will 
give  the  best  results.  To  determine  velocity  information  from 
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the  variance,  skewness,  or  kurtosis  of  the  translation  scans, 
several  scans  at  3,  4,  or  5  angles  respectively,  are  necessary.  A 
better  method  is  introduced  in  the  next  chapter  which  studies 
the  scattering  phenomenon  as  a  function  of  angle. 


APPLICATION  OF  THE  SPECTRAL  THEORY  TO  THE 
CASE  OF  VARYING  ANGLE  AND  FIXED  SCAN 

POSITIONS 


1. 1 


The  second  experiment  applies  the  spectral  theory  to  the 
case  of  rotating  the  crossed  sound  beams  at  a  single  scan 
position.  From  each  360*  scan,  turbulent  velocity  correlations 
can  be  determined  similarly  to  the  method  outlined  for 
translational  scans.  But  to  solve  a  system  of  86  equations  (one 
for  each  angle  scanned)  in  closed  form  is  inefficient.  Instead 
two  numerical  methods  are  proposed  to  derive  the  velocity 
information.  One  is  based  upon  Fourier  analysis,  and  the  other 
is  based  upon  a  linear  least  squares  fitting  algorithm. 

The  Fourier  method  re-interprets  the  Doppler  shift 
equations  for  a2,  S,  K  and  <  f  >  as  finite  Fourier  series. 

An  example  of  this  method  (using  the  variance  from  Eq.  11)  is 
outlined  below: 


a2(0.)  =  KX2<VX2>  +  2KxKy<VxVy>  +  Ky2<Vy2> . 
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The  quantities  Kx2,  KxKy,  Ky2,  and  a2  are  all  functions  of  8.. 
However,  Kx2,  KxKy,  and  Ky2  can  be  written  as  a  finite  Fourier 
series  in  8*  without  too  much  difficulty.  Start  with  the  Kx2 
expression  below: 

Kx2  =  ko+2(l  -  ^cos  0,  -  ysin  0.)2 


This  expression  can  be  expanded  and  regrouped  using 
trigonometric  identities  to  yield: 


K*2  =  k0+2[ 


v4  2. 

y  ^sin  28.  ]  . 


VI  cos  8.  + 


2  Icos  28.  -  n  sin  8.  + 


(21) 


Similarly,  Fourier  expressions  can  be  found  for  KxKy  and 
Ky2  below: 

KxKy  =  ko*2[  frg  +  T)  .  VlEcos  e.  +  (y-^cos  26.  -  fsin  6.  + 

'  4  z  4  z 

sin  20,],  (22) 

.  (l+Y2)  (y2-1)  Y 

Ky2  ~ko+2[- - L  +  \ - Zcos  20*  +-sin  20*].  (23) 

4  4  2 

Thus  o2  can  be  rewritten  as 


or2  =  ko+2[ 


4  "  T  <vx2>  +  2^y^  +  y^  <VxVy>+  <Vy2>] 


+  ko+2[-  VI  <VX2>-  2y  ^<VxVy>]  cos  0. 


+  fco+2[ 


1  £ 
v4  "  2  2 


<vx2>  +  2 <VXV '^<V  2>|cos28. 
4  J  4  y 

+  ko+2[  -  ZY<VX2>  -  2^  <VxVy>]sin  0. 

+  k0+2[Y  ^  <VX2>  +  2  V-+  <VxVy>+ 1< Vy2>]sin  20.  ,  (24) 


which  is  now  a  Fourier  expansion  of  a2-  A  Fast  Fourier 
Transform,  FFT,  can  be  performed  upon  the  a2  data.  Although 
the  a2  experimental  data  does  not  contain  2n  points,  which  is  a 
prerequisite  for  performing  FFTs,  the  data  repeats  every  full 
360*  revolution  of  the  transducers.  Therefore,  one  can  expand 
the  data  to  2n  points  by  adding  repetitions.  Here,  n  denotes  an 
integer  value.  From  the  FFT,  one  will  find  the  coefficients  (a,  b, 
c,  d,  and  e)  of  a2. 


a2  =  a  +  b  cos  0.  +  c  cos  20.  +  d  sin  0.  +  e  sin  20. ,  (25) 

which  must  be  equivalent  to  the  coefficients  of  the  Kx  and  Ky 
terms.  This  is  true  mathematically  using  the  orthogonal 


property  of  sines  and  cosines.  This  sets  up  the  following 
system  of  equations: 
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a 

=  ko.2l 

(l.f 

4  2 

V 

j<Vx2>  +  2^—  *  Y)<vxv>>+  +  ^^Vy2^ 

(26) 

b 

=  ko+2[-  VI  <VX2>-  27  ^<VxVy>] 

(27) 

c 

<N 

+ 

O 

M 

II 

1  ±\ 
4  2 

<VX2>  +  2^— <VxVy>  +  ^  '  ^<Vy2>] 

(28) 

d 

=  ko*2l 

-  Tt  <VX2>  -  21?-  <VxVy>] 

(29) 

e 

=  ko.2n  f  <VX2>  +  2 <vxvr>+  J<Vy2>] . 

(30) 

From  this  system,  one  can  solve  to  get  <VX2>,  VxVy>,  and  <Vy2>. 

Similar  finite  Fourier  sums  have  been  calculated  for  <CDd>, 
<C0d  -  <C0d>)3>,  and  <(C0d  -  <C0d>)4>.  These  sums  can  be 
interpreted  similarly  to  the  variance  data  to  find  V0x,  Voy, 

<V„3>,  <Vx2Vy>,  <V;Vyb..  <vyV, 

<VX4>,  <VX3V)>,  <Vx2Vy2>,  <VxVy3>,  and  <Vy->. 

The  other  method  of  finding  the  velocity  correlations  uses 
an  algorithm  known  as  the  method  of  least  squares2.  This 
procedure  involves  fitting  a  function  F  to  the  experimental 
data.  Adjustable  parameters  of  the  function  F  are  determined 
by  minimizing  an  error  function,  E,  as  described  in  the 
following  example: 

In  Equation  11,  the  variance,  a2,  is  defined  as: 


a2(6.)  =  KX2<VX2>  +  2KxKy<VxVy>  +  Ky2<Vy2>. 

Let  f  (0)  =  Kx2(0) 

g(0)  -  2KxKy(0) 
h(0)  =  Ky2(0) 
a  =<VX2> 
b  =  <VxVv> 

C  = 

and  F(a,  b,  c,  0)  =  af(0)  +  bg(0)  +  ch(0) . 
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The  error  function  is  E  =  X  {o2-F}2.  (31) 

all  data  pis 

Now  the  adjustable  parameters  in  this  problem  are  the  velocity 
correlations  a,  b,  and  c.  If  they  are  chosen  appropriately,  the 
error  function  will  be  a  minimum.  This  means  that  the 
calculated  velocity  correlations  are  the  values  the  theory 
predicts  for  that  data. 

To  find  the  minimum  errors,  find  the  minimum  of  E  with 


respect  to  a,  b,  and  c. 

The  minimum  for  a  is 

dE 

&  =  ( X  die2  -  F)2)  /  3a 

=  pto2  -  F)(-f(0)).  (32) 

Set  derivative  equal  to  zero  to  find  extrema 
I  f(8)a2  =  £f(6)  F 

=  a£  f2(9)  +  b£  f(9)g(8)  +  c£  f(9)h(0).  (33) 

Similarly  the  minima  of  b  and  c  yield  the  relations 

X  8(8) o2  ■  *X  f(0)g<9)  +  b£  g2(0)  +  c£  g(9)h(0)  ,  (34) 

X  h(0)o2  =a£f(0)h(0)  +  b£  g(0)h(9)  +  c£  h2(0) .  (35) 

Eqs.  33,  34,  and  35  yield  the  following  system  of  equations: 

=  a  Bff  +  b  Bfg  +  c  Bfh  (36) 

Cgo  =  a  Bfg  +bBgg  +cBgh  (37) 

=  a  Bn,  +  b  Bgh  +  c  Bhh  ,  (38) 

where 

Cfo  =  X  K,2(0)o2 
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C«°  =X  2^(6)^ 

Cte  =  X  Ky2(0)c2 

Bff  -  X  *<x2(e)Kx2(fl) 

Bfg  =  X  2Kx2(e)KxKy(0) 

Bfh  =  X  K,2(e)Ky2(0) 

Bfg  =  X  4KxKy(e)KxK/e) 

Bgh  =  X  2KxKy(6)Ky2(6) 

Bhh  =  X  Ky2(8)Ky2(6) .  (39) 

This  system  of  equations  can  be  solved  by  Gauss-Jordan 
row  reduction  which  is  easily  handled  by  a  computer.  The 
computer  routine  to  perform  Gauss-Jordan  reduction  has  been 
written  for  the  Macintosh  II  by  the  author.  This  program  is 
listed  in  Appendix  C. 


REYNOLDS  NUMBER  SIMILARITY  IN  TURBULENT 

FLOW 


In  Townsend’s  classic  text,  The  Structure  of  Turbulent 
Shear  Flow.9  on  page  89  he  points  out 

"that  a  turbulent  flow,  suitably  defined  by 
boundary  conditions,  has  a  main  structure 
that  is  independent  of  the  fluid  viscosity, 
provided  that  it  is  not  too  large." 


He  adds 


"...  geometrically  similar  flows  are  similar  at 
all  sufficiently  high  Reynolds  numbers..." 

Finally,  Townsend  states  the  principle  of  Reynolds  number 
similarity. 


"In  a  fully  turbulent  flow,  there  exists  a 
region  including  almost  all  the  flow,  over 
which  the  direct  action  of  viscosity  on  the 
mean  flow  is  negligible,  ..." 

Commenting  that  the  Reynolds  stresses  are  large  compared 
with  the  mean  viscous  stresses,  he  states  that  the  energy- 
containing  components  of  the  turbulence  are  determined  by 
the  boundary  conditions  of  the  flow  alone  and  are  independent 
of  the  fluid  viscosity. 

We  have  made  use  of  this  fact  in  our  project. 
Wygnanski12  studied  turbulent  flow  in  air  at  Reynolds 
numbers  on  the  order  of  Re  =  105.  Thus,  his  data  can  be 
compared  to  our  own.  From  fits  and  trends  in  our  data,  we  will 
be  able  to  decide  how  well  we  predict  turbulent  flows. 
Wygnanski's  results  are  considered  the  standard  of 
measurement  for  turbulent  jet  flow  in  air.  Thus,  if  we  can 
show  results  similar  to  his,  we  will  have  succeeded  in  proving 
our  method  actually  measures  turbulent  velocities.  However, 
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we  should  encounter  differences  in  the  higher  order  correlation 
functions  because  the  physical  make-up  of  the  jet  used  by 
Wygnanski  and  the  one  used  in  the  scattering  experiment  is 
different.  Wygnanski  uses  a  simple  nozzle  (of  diameter  1.04 
inches)  to  generate  turbulent  shear  flow  in  a  large  open  space. 

Our  submerged  water  jet  has  a  large  circular  disk  (of  diameter 
equal  to  6  inches)  encompassing  the  1/4  inch  nozzle  and  is 
placed  in  a  one  cubic  meter  tank.  Unfortunately  the  tank  is  not 
large  enough  to  avoid  interfering  with  the  entrainment  of  the 
flow  on  the  skirts  of  the  jet.  This  setup  was  unavoidable  due  to 
size  constraints  of  our  laboratory  and  tank  availability. 
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EXPERIMENTAL  RESULTS 

A.  Description  of  the  Spectral  Analysis  Measurements 

The  results  of  this  experiment  present  the  calculated 
velocity  distributions  across  the  width  of  the  turbulent  water 
jet.  The  process  by  which  measurements  of  the  nonlinearly 
scattered  sum  frequency  intensity  spectra  are  used  to  obtain 
the  turbulent  velocity  correlations  involves  several 
computational  steps  that  are  explained  next.  Translational  and 
rotational  data  analysis  invole  similar  discussions.  Therefore  a 
single  general  discussion  is  presented  for  both  cases. 

It  is  pointed  out  that  spectra  (measured  on  the  Tektronix 
7L5  sweep  spectrum  analyzer)  are  a  detection  of  the  sum 
frequency  tone  that  changes  its  frequency  as  a  function  of  time. 
Therefore,  one  must  average  the  results  of  many  spectral 
sweeps  in  order  to  obtain  a  true  time  average  spectrum.  The 
IQ400  digital  oscilloscope  accomplishes  this  task  by  signal 
averaging  the  output  voltage  from  20  sweeps  versus  time.  The 
output  voltage  is  proportional  to  the  rms  pressure  amplitude  in 
the  spectrum.  The  sweep  time  base  is  proportional  to  the 
frequency  in  the  spectrum.  The  average  rms  spectrum  is  then 
saved  on  a  3.5  inch  floppy  disk.  From  voltage  and  time  base 
calibrations  the  average  rms  spectrum  is  converted  to  an 
intensity  spectrum  as  a  function  of  frequency.  The  averaged 
spectra  for  several  angles  are  shown  in  Figs.  19a  and  19b.  Each 
spectrum  was  obtained  with  the  overlap  region  fixed  upon  the 
jet  axis. 

The  translational  scattering  experiments  (at  fixed  angle) 
involved  scanning  the  overlap  region  at  40  scan  positions 
across  the  jet  over  the  range  x=  ±  2  inches.  The  angular 
scattering  experiments  are  performed  at  80  individual  angles 
at  each  scan  position.  Angular  scattering  was  performed  at  1 1 
different  scan  positions  between  x=0  and  x=  2.0  inches.  The 
entire  crossed  beam  scattering  experiment  involved  880  time 
averaged  spectra. 


B.  Angular  Scattering  Results 


The  graphs  of  the  total  intensity  I+(8,)  ,  Doppler  shift  fd  = 
<  f  >  -  fo+  ,  variance  a2,  skewness  S,  and  kurtosis  K  are 
presented  as  a  function  of  angle  0.  for  the  interaction  region 
located  on  the  jet  axis  (Figs.  20  -24).  The  calculations  for  the 
spectral  moments  for  the  average  spectrum  are  obtained  by  a 
computer  analysis.  An  Apple  He  Basic  program  called 
"Analyzer"  is  used  to  analyze  the  translational  scans  and  a 
Turbo  Pascal  program  called  "Mac"  is  used  for  the  angular 
scattering  trials.  These  programs  converted  the  rms  voltage 
spectra  to  intensity  spectra  and  then  found  the  spectral 
moments.  These  programs  are  listed  in  Appendix  C. 

In  order  to  obtain  the  required  velocity  correlations,  the 
values  from  the  spectral  moments  <  f  >,  <(f  -  <  f  >)2>, 

<(f  -  <  f  >)3>,  and  <(f  -  <  f  >)4>  are  numerically  curve  fit  using  the 
linear  least  squares  fit  algorithm  that  was  outlined  on  pages  32 
-  33.  The  results  for  the  Doppler  shift,  fd,  are  shown  in  Fig.  21. 
The  results  for  the  variance,  a2,  (or  second  moment)  along  the 
third  and  fourth  moments  are  shown  in  Figs.  22,  25  and  26, 
respectively. 

Results  for  determining  the  mean  velocity  components. 
Vox  and  V0y,  involve  measurements  of  the  Doppler  shift,  fd  = 
(<  f  >  -  fo+)  across  the  jet.  These  velocities  (shown  in  Figs.  27 
and  28)  are  plotted  versus  the  normalized  scan  position  "H, 
which  is  given  by  x  /  y  .  Here,  the  scan  position  is  located  at 
the  distance  y  =  33.9  d  from  the  nozzle  exit,  where  d  is  the 
nozzle  diameter. 

Pitot  tube  measurements  (of  the  axial  mean  flow  velocity 
profile  across  the  jet)  obtained  in  the  laboratory  are  used  to 
determine  the  accuracy  of  the  nonlinear  scattering  method.  A 
comparison  of  the  pitot  tube  measurements  of  Voy  across  the 
jet  with  the  rotational  scattering  experiment's  determination  of 
Voy  is  shown  in  Fig.  27.  On  axis,  where  “H  =  0,  the  pitot  tube 
maximum  axial  velocity,  Um  ,  is  1.187  m/s.  For  comparison,  the 
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the  scattering  experiment  yielded  a  value  of  Um  =  1.22  ±  0.01 
m/s.  These  profiles  of  the  axial  mean  flow  are  in  agreement 
with  each  other  and  with  Wygnanski’s12  published  data  (not 
shown)  for  a  similar  jet  flow  involving  a  turbulent  air  jet. 

The  radial  mean  velocity  distribution  Vox  is  shown  in  Fig. 
28.  This  velocity  profile  is  theoretically  obtained  by  inserting 
the  axial  mean  velocity  Voy  into  the  equation  of  continuity  for 
mean  flow.  Notice  that  Vox  is  on  the  order  of  10'2  times  less 
than  Voy,  making  it  difficult  to  predict  Vox  front  experimental 
scattering  data.  Pitot  tube  measurement  of  the  radial  flow  are 
not  possible. 

The  second  set  of  graphs  are  related  to  the  rms  velocity 
distributions.  Since  translational  scattering  data  and  pitot  tube 
measurements  are  unable  to  obtain  the  higher  order 
correlations,  Wygnanski's  data12  is  used  to  make  comparisons 
with  the  rotational  measurements  of  the  rms  velocities.  The 
predictions  for  the  second  order  correlation  <VX2>  (Fig.  29) 
matched  well  with  Wygnanski’s  data.  Wygnanski's  jet, 
however,  does  not  have  the  6  inch  diameter  plate  surrounding 
its  nozzle.  This  might  account  for  the  reason  our  radial 
variations  are  broader.  The  other  two  second  order  velocity 
correlations  <VxVy>  and  <Vy2>  (shown  in  Figs.  30  and  31), 
followed  the  same  general  trend  as  Wygnanski,  but  the  results 
shown  here  indicate  a  subtle  difference  between  his  jet's 
construction  and  ours. 

The  third  order  velocity  correlation  coefficients  are 
presented  in  Figs.  32-35.  These  values,  determined  from  the 
rotational  runs,  are  compared  with  some  of  Wygnanski’s  data. 

In  Fig.  33,  it  is  observed  that  both  sets  of  plots  for  <VX3>  agree 
in  magnitude  and  trend.  This  is  also  true  for  the  results  of 
<Vx2Vy>  that  are  shown  in  Fig.  32.  However,  a  discrepancy 
arises  with  the  <VxVy2>  profile  of  Fig.  32  which  indicates  that 
the  differences  in  his  jet  and  ours  is  noticeable  at  higher  order 
velocity  correlations. 

The  fourth  order  velocity  correlations  are  presented  in 
Figs.  36  and  37.  These  correlations  are  a  triumph  of  the  project 


because,  to  the  author's  knowledge,  no  one  else  to  date  has 
been  able  to  measure  them  reliably.  Since  there  are  no 
published  fourth  order  relations  to  compare  with  our  jet,  we 
cannot  make  any  comparisons.  Structure  can  be  observed  in 
the  functions.  This  indicates  that  the  apparatus  is  sensitive  to 
fourth  order  fluctuations. 

The  last  graphs  shown  in  Figs.  38  and  39  depict  the 
variations  in  turbulent  intensity  across  the  jet.  One  defines  the 
turbulent  intensity  to  be  the  ratio  of  the  rms  fluctuating  motion 
in  a  specified  direction  to  the  mean  flow  measured  at  that 
point.  The  intensity  profiles  show  the  dramatic  increase  in 
turbulent  rms  velocity  off  axis.  These  results  indicate  that  the 
flow  becomes  increasingly  random  off  axis. 

C.  Translational  Scanning  Results 

Nonlinear  scattering  experiments  that  are  performed  at  a 
fixed  angle,  while  the  interaction  region  scans  across  the  width 
of  the  jet,  are  called  translation  scans.  For  the  experiments 
involving  translation  scans  at  fixed  angle  (forward  scattering  at 
0*  and  back  scattering  at  180*).  Predictions  for  the  radial  and 
axial  mean  velocity  components,  Vox  and  Voy  can  be  obtained 
by  measuring  the  average  frequency  <  f  >  vs  scan  position  (and 
subsequently  the  Doppler  shift  fd  =  <  f  >  -  f0+)  and  directly 
inserting  these  results  into  the  closed  form  solutions  given  by 
Eqs  19  and  20.  The  Doppler  shift  as  a  function  of  scan  position 
"H  (for  0.  =  0  and  0*  =  180  )  are  shown  in  Fig  40  and  41, 
respectively.  These  Doppler  results  are  used  in  Eq.  19  to 
predict  the  radial  mean  flow  velocity  profile,  Vox(*l),  that  is 
shown  in  Fig.  42. 
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FIG.  19b.  Sum  frequency  spectra  for  angles  -117 
-30  degrees  and  30  to  117  degrees 
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FIG.  21.  Regression  plot  of  the  average  Doppler  shift 
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FIG.  22.  Regression  plot  of  the  variance 
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FIG.  23.  Skewness  vs.  theta  star 
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24.  Kurtosis  vs.  theta  star 
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FIG.  26.  Regression  plot  of  the  fourth  moment 
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FIG.  29.  Comparison  of  Wygnanski’s  data  to  NL  X-Beam 
data  (rotational)  for  the  radial  rms  velocity  across  the  jet 
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FIG.  32.  Triple  Velocity  Correlations 
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FIG.  33.  Comparison  of  Wygnanski's  data  to  NL  X-Beam 
data  (rotational)  for  the  third  order  radial  velocity 
component 
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FIG.  35.  Triple  velocity  correlation  coefficient  (axial)  across 
the  jet 
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FIG.  36.  Quadruple  velocity  correlations  across  the  jet 
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FIG.  39.  Distribution  of  the  axial  turbulent  intensity 
across  the  jet 
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DISCUSSION  OF  POSSIBLE  ERRORS 

We  are  able  to  map  out  the  turbulent  velocities  for  a 
submerged  water  jet  with  the  nonlinear  scattering  apparatus. 

This  statement  is  founded  upon  the  agreement  between  the 
results  obtained  by  our  analysis,  general  comparisons  with  the 
measurements  obtained  by  Wygnanski,  and  by  comparisons 
with  the  pitot  tube  results.  The  disagreement  in  the  detailed 
nature  of  the  profiles  of  the  higher  order  correlations  is 
attributed  to  the  differences  in  the  turbulent  jets.  Wygnanski 
uses  2  bare  nozzle  in  a  large  room.  Our  jet  nozzle  is  baffled  by 
a  6  inch  diameter  flat  plate  and  ejects  into  a  cubic  meter  tank 
of  water.  The  boundaries  of  the  sides  of  the  tank  and  opposite 
wall,  in  addition  to  the  entrained  flow  having  to  go  around  the 
baffled  orifice,  are  suspected  to  have  interfered  with  a  more 
simple  entrainment  of  fluid  into  the  jet.  These  differences  are 
known  to  have  small  effects  upon  the  structure  of  the  mean 
flow  profiles.  The  second  order  correlations  were  not  expected 
to  compare  well  with  Wygnanski’s  jet  -  although  they  did.  At 
the  higher  orders,  the  correlations  do  generally  compare  well, 
but  not  in  detail.  It  is  believed  that  this  scattering  experiment 
can  predict  mean  flows  to  «  1.5  %,  second  order  correlations  to 
*  5  %,  third  order  correlations  to  »  10  %,  and  fourth  order 
correlations  to  *  25  %,  if  the  turbulence  is  not  too  weak. 

There  are  some  large  errors  in  the  tails  of  our  curves  ( "n  > 
0.15)  that  were  caused  when  radio  frequency  interference 
from  the  Severn  Naval  Station's  radio  transmissions  was 
detected  in  the  sum  frequency  spectra.  The  interference  forced 
us  to  delete  some  of  the  spectra  that  were  obtained  in  the  set 
of  80  scattering  angles  per  scan  position.  This  in  turn  reduced 
the  number  of  data  points  that  were  used  to  fit  the  curves. 

Therefore,  our  curve  fitting  routines  on  the  skirts  of  the  jet 
(having  fewer  data  points)  resulted  in  larger  errors  in  the 
predictions  of  the  turbulent  velocity  correlations 

It  is  noticed  that  near  the  edges  of  the  jet,  errors  are 
generally  larger  for  the  axial  velocity  measurements  than  the 


radial  ones.  This  may  be  explained  in  terms  of  the  optimum 
scanning  angles.  At  the  optimum  angles,  correlations  between 
the  spectral  moments  and  the  turbulent  velocities  leads  to  a 
least  error.  Angular  scattering  for  0.  between  60*  and  120’,  or 
between  240*  and  300*  are  not  possible  as  the  transducer  units 
are  translated  off  axis.  Here,  the  transducer  units  are 
physically  moved  into  the  turbulent  jet.  Unfortunately,  the 
maximum  values  of  Ky  occur  in  this  range.  Doppler  shift 
contributions  from  the  factor  KyVy  need  to  be  large  if  the 
velocity  component  Vy  to  be  measured  is  very  small. 

Therefore,  reducing  the  scans  in  this  range  degrades  the  ability 
of  the  apparatus  to  determine  the  axial  flow.  The  angles  from 
0.  =  -30*  to  30*,  and  from  150*  to  210*  are  the  extreme  forward 
and  back  scattered  angles.  In  translating  the  interaction  region 
across  the  jet  at  these  angles,  the  transducer  units  do  not 
interfere  with  the  turbulent  jet.  Therefore,  a  full  range  of 
values  for  Kx  is  covered  at  each  radial  position  producing 
Doppler  shift  contributions  from  the  factor  KXVX  that  are 
significant.  Moreover,  this  extreme  forward  or  back  scattering 
range  corresponds  to  a  maximum  for  Kx.  Therefore,  Doppler 
shift  predictions  from  the  factor  KXVX  lead  to  very  small  errors 
in  predicting  Vx. 


CONCLUSIONS 
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Nonlinear  crossed  beam  scattering  in  the  presence  of 
turbulence  produces  a  sum  frequency  component  that  radiates 
outside  the  interaction  region.  Measurements  of  the  Doppler 
shift,  spectral  broadening,  skewness,  and  kurtosis  (of  the 
scattered  intensity  spectrum)  can  be  used  to  predict  mean  flow 
velocities  as  well  as  second,  third,  and  fourth  order  turbulent 
velocity  correlation  coefficients.  The  performance  of  the 
crossed  beam  experiment  in  measuring  turbulent  flow  is  very 
sensitive  to  the  ratio  ^  =  ko-  /  ko+.  If  “V  «  1,  then  scattered 
measurements  in  the  forward  (0,  =  0‘)  or  back  (0.=  180') 
scattered  directions  will  be  highly  correlated  with  the 
component  of  velocity  along  the  direction  of  the  bisecting  line 
between  the  sending  transducer  axes  and  the  point  of  overlap. 

This  nonlinear  scattering  phenomenon  is  highly  sensitive 
to  even  the  smallest  vorticity  disturbances  that  were  created  in 
the  water  tank.  This  diagnostic  tool  is  a  viable  alternative  to 
conventional  turbulent  flow  measuring  devices  that  use  hot 
film  probes  or  laser  Doppler  techniques. 
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A.  Single  beam  scattering  of  sound  by  turbulence 

The  scattering  of  sound  by  turbulence  predicts,  in  theory, 
that  the  scattered  sound  waves  will  be  Doppler  shifted  by  an 
amount  equal  to  K  •  V,  where  K  =  ks  -  k;  is  the  wavevector  of 
the  scattered  wave,  ks,  minus  the  wavevector  of  the  incident 
wave,  kj.  Define  the  velocity  vector  of  a  fluid  particle  to  be  V. 
One  can  derive  the  Doppler  shift  expression  by  studying  the 
relative  phase  of  an  incident  acoustic  wavefront  that  has  been 
scattered  by  a  moving  fluid  particle  or  turbulent  eddy.  Figure 
43  illustrates  how  a  turbulent  eddy  can  scatter  sound. 


FIG.  43.  Scattering  of  an  incident  sound  wave  by  a  turbulent 

eddy 

If  the  wavefront  is  deflected  through  an  angle  6  by  the 
action  of  the  moving  fluid  particle,  there  is  a  change  in 
direction  of  the  wavefront.  The  change  in  the  wave’s  path 
causes  one  portion  of  the  wave  front  to  travel  a  longer  distance 
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(ab+bc)  than  the  other.  After  the  scattering,  the  waves 
propagate  along  the  new  path.  However,  because  of  this  path 
difference,  the  scattered  wave  will  exhibit  a  phase  shift  (with 
respect  to  the  incident  wave)  that  is  given  by: 

A<j>  =  kj  ab  +  ks  be. 

For  turbulent  fluid  particles,  the  particle  velocity  V,  is  always 
changing  in  time.  This  causes  the  relative  paths  ab  and  be  to 
vary  with  time.  These  path  variations  change  the  phase  shift, 
A<(>  with  respect  to  time.  One  can  relate  the  time  dependent 
phase  shifts  in  the  scattered  wave  to  the  Doppler  shifts,  a>d, 
using  the  expression  below. 

A<J> 

—  =  cod  =  -kjv  sin  x  -  ksv  sin  x  .  (Al) 

At 

A<t> 

Here,  ~  is  the  time  rate  of  change  in  phase.  The  negative  signs 

appear  in  this  expression  because  the  particle  velocity 
components  move  away  from  the  source  and  away  from  the 
receiver,  respectively.  This  results  in  a  net  down  Doppler  shift. 
If  kj  =  ks  then  ©d  becomes 

tod  =  -  2ks  v  sin  (0/2),  where  0  =  2x.  (A2) 

Now  we  replace  v  in  the  above  model  with  V,  the  vector 
turbulent  particle  velocity.  The  magnitude  of  ab  is  described 
by 


ab  =  k;  V  cos  a  =  k;  *  V.  (A3) 

The  value  is  negative  because  the  direction  of  V  along  k,  points 
away  from  the  transducer.  Similarly, 

be  =  ks  Vcos  a  =  ks  •  V  .  (A4) 
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This  value  is  positive  since  be  is  towards  the  receiver.  The  net 

=  -  (k*  •  V)  +  (iT  .  V) 

=  <£-?)•  V 

=  K  •  V  .  (A5) 

B.  Nonlinear  crossed  beam  scattering  of  sound  by  turbulence 

It  is  known  from  Korman  and  Beyer’s  nonlinear 
scattering  theory5  that  the  nonlinearly  scattered  sum 
frequency  component  also  undergoes  a  Doppler  shift.  Using 
this  model,  one  can  relate  ks  and  kj  to  nonlinear  crossed  beam 
parameters.  Let  the  nonlinearly  scattered  sum  frequency 
wavevector  ks  be  equal  to 

ks  =  ko+  n  =  (koi  +  ko2)  n  ,  (A6) 

where  n  is  a  unit  vector  that  points  from  the  interaction  region 
to  the  receiver.  Define  the  incident  wavevector,  as  the  linear 
combination  of  the  incident  wavevectors  from  each  primary 
wave.  Then,  kj  becomes 


result  is 


A0 

At 


=  ©d 


K  ~  koi  n0i  +  k02  n02  t  (A7) 

which  are  the  two  emitted  waves  from  the  transmitting 
transducers  #1  and  #2.  The  difference  between  the  scattered 
and  incident  wave  vectors,  Ak  ,  can  be  expressed  by 

*■"  * 

k5 '  k;  =  Uc01  k01  n01]  +  [k02  n  -  k02  n02]  , 

=  K,  +  K2 

=  K-  (A8) 
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Here,  one  has  defined  K}  and  K2  as  individual  scattering  vectors 
of  the  nonlinear  scattering  experiment.  Further,  define  K+  =  Kj 
+  K2  to  be  the  resultant  scattered  wave  number  for  the 
nonlinear  scattering  theory.  Thus  we  have  our  relation  for  the 
Doppler  shift,  o)<j  =  K+  •  V  which  is  based  upon  the  scattering  of 
longitudinal  sound  waves  from  transverse  (or  rotational)  fluid 
particle  flow. 

C.  Doppler  quantum  theory 

The  nonlinear  scattering  phenomena  involves  four 
principle  parts:  two  incident  focused  sound  waves,  a  scattered 
wave,  and  a  turbulent  eddy  that  scatters  the  incident  wave.  To 
examine  the  scattering  effects  in  the  interaction  region,  first 
study  a  model  for  conventional  single  beam  scattering.  Later, 
one  can  expand  the  model  to  include  both  incident  beams  and 
study  the  nonlinear  scattering  from  crossed  beams. 

1.  Conventional  single  beam  scattering  theory 

In  conventional  single  beam  scattering,  an  incident  sound 
wave  is  scattered  by  the  action  of  a  turbulent  eddy.  The 
resulting  momentum  change  is  shown  in  Figure  44. 


k 

1 


q 


FIG.  44.  Momentum  change  of  incident  beam  due  to  scattering 
To  find  a  relationship  between  the  incident  and  scattered 
waves,  use  the  principles  of  momentum  and  energy 
conservation.  The  momentum  and  energy  of  an  excitation  (or 
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wave  packet)  may  be  expressed  by  hk,  ho),  respectively,  where  h 
is  Planck's  constant  divided  by  2x,  k  is  a  wave  vector  of  a 
phonon  excitation,  and  to  is  the  angular  frequency  of  the  wave 
excitation.  Thus,  the  momenta  of  the  incident  and  scattered 
phonons  become  hkj  and  hks  ,  respectively.  The  energies  of  the 
phonon  excitations  are  hcoj  and  hcos.  For  the  scattering  model, 
let  q  denote  the  wave  vector  excitation  of  the  turbulent  eddy 
and  o)t  the  corresponding  angular  frequency  of  the  turbulent 
eddy  excitation.  Thus,  its  momenta  and  energy  are  hq  and  hcot. 

The  conservation  principles  are 

kj  +q  =  ks  conservation  of  momentum  (A9) 

©i+cot  =  <os  conservation  of  energy  .  (A10) 

From  momentum  and  energy  conservation,  one  can  solve  for 
the  scattered  wave  number  q  and  the  turbulent  excitation 
frequency: 

q  (All) 

cot  =(cos-<Bi).  (A  12) 

Define  q  to  be  K  ,  the  scattered  wave  vector,  and  a>t  to  be  to  a  , 

the  Doppler  shift.  Thus  the  excitation  frequency  of  the 
turbulent  eddy  is  equal  to  the  Doppler  shift  in  the  incident 
wave.  For  sound  waves,  the  wave  number  is  directly  related  to 
its  excitation  frequency. 

<*>i  =  ki  c  (A  1 3) 

t»s  =ksc,  (A  14) 

where  c  is  the  speed  of  sound  in  water.  For  a  turbulent  eddy, 
the  excitation  frequency  is  known  to  be 

=q  •  V,  (A  1 5) 

where  V  is  the  velocity  of  the  turbulent  eddy. 


Make  the  following  convenient  substitutions: 
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kj  =  ki  iij  ,  (A16) 

where  kj  is  the  wave  number  and  ni  is  a  unit  vector  in  the 
direction  of  propagation  of  the  incident  wave.  Similar 
substitutions  are  used  for  the  scattered  wave: 

K  -  ks  ns  .  (A17) 

From  the  equation  of  momentum  conservation 

q  =  (ks  ns  -  ki  nO  ,  (A18) 

the  turbulent  excitation  frequency  is  now  given  by 

©t  =q  •  V 

=  (ksns  -  kj  n i)  •  V  .  (A19) 

The  Doppler  shift  is  the  excitation  frequency  wt  or  ©s  -  ©i  .  It 
can  be  expressed  by 

©s-wi  =  (©s  ns  -  ©j  nj)  •  V/c  .  (A20) 

The  Doppler  shift  is  now  cast  in  a  form  that  can  be  used  to 
interpret  the  experimental  geometry.  Now,  one  can  solve  for 

©s: 


©s  =--  ©j[l  -  nj  •  (V/c)]  /  [1  -  ns  •  (V/c)]  .  (A21 ) 

In  the  experiments  of  interest,  the  magnitude  of  the  flow 
velocity  I  V  I  is  much  less  than  the  speed  of  sound,  c. 

Since  I  V  l/c  «  1 , 

©s  =©i[l  -  nj  •  (V/c)][  1  +  ns  •  (V/c)]  .  ( A22) 


If  one  drops  the  second  order  term  in  A22,  the  relation 
between  the  scattered  wave's  angular  frequency  and  the 
incident  wave's  angular  frequency  becomes 

G>s  ~t»>j[l  +  (ns-nj)«  (V/c)]  .  (A23) 

These  relations  predict  that  the  scattered  waves  will 
gain/lose  energy  {depending  upon  the  sign  of  (ns  -  nj)  4  (V/c)} 
from  the  excitations  which  have  been  modeled  as  fluctuations 
of  the  turbulent  eddies. 


Up  Doppler  shift 


FIG.  45.  Scattering  geometries  that  produce  up  Doppler  or  down 

Doppler  shifts 

In  Figure  45,  the  scattering  diagram  on  the  left  describes  a 
scattering  geometry  in  which  the  (ns  -  nj)  •  (V/c)  term  is 
positive,  representing  an  up  Doppler  shift.  Conversely,  the 
diagram  on  the  right  shows  a  scattering  geometry  in  which  the 
(ns  -  nj)  •  (V/c)  term  is  negative,  corresponding  to  a  down 
Doppler  shift. 
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2.  Nonlinear  crossed  beam  theory 

In  nonlinear  crossed  beam  scattering,  two  incident  sound 
beams  cross  in  the  presence  of  a  turbulent  eddy.  The 
interaction  creates  a  scattered  sum  frequency  component 
whose  frequency  is  the  algebraic  sum  of  the  frequencies  of  the 
incident  beams.  The  resulting  momentum  change  for  this 
interaction  is  shown  in  Figure  46. 


FIG.  46.  Momentum  changes  of  incident  beams  due  to  nonlinear 

scattering 

In  the  nonlinear  interaction,  there  are  two  distinct 
incident  waves,  each  with  its  own  momentum  and  energy.  As  a 
result,  the  conservation  equations  become 

k,  +  k2  +  q  =  ks,  (A24) 

C0i  +  0)2+G)t  =  CDS.  (A25) 

One  can  solve  equations  A24  and  A25  to  find  a  relationship 
between  q  and  tot: 


q 


=  ks  -  (k,  +  k2) 


(A26) 


©t  =  ©s  -  (0)1  +  ©2)  • 


(A27) 


These  relations  become  more  useful  by  making  the  following 
definitions: 

<n+  =  ©1  +  ©2 

kj  =  koi  not 

I4  =  ko2  no2 

ks  =  (koi  +  ko2)n  .  (A28) 

Here,  define  two  wave  numbers 

koi  =  ®oi  /  c  , 
ko2  =  ©02  /  c  , 

and  their  corresponding  incident  unit  vector  directions  to  be 
noi  and  n<)2  .  respectively.  Let  n  be  a  unit  vector  in  the 
direction  of  the  scattered  wave.  From  these  definitions  q  and  ©t 
become 


q  =  (koi  +  ko2  )n  -  (koi  "01  +  ^02  "02)  (A29) 

©t  =©s-©+,  (A30) 

which  can  be  cast  in  the  form 

q  =K  (A31> 

©  =  ©s  -  ©+  (A32) 


by  using  the  following  definitions 

koi  (n-noi) 
ko2  (n-no2) 
icj  +  K^  . 


K,  = 

^2  = 

K  = 


(A33) 
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Thus,  the  resulting  Doppler  shifts  act  as  if  the  superposition  of 
two  simultaneous  conventional  single  beam  scatterings  had 
occurred. 

From  the  conservation  equations  for  nonlinear  scattering, 
one  can  prove  that  two  collinear  beams  can  produce  a  radiated 
sum  frequency  component  whereas  two  crossed  beams  cannot 
unless  turbulence  is  present  in  the  interaction. 

In  the  absence  of  turbulence,  q  and  tot  become  zero, 
leaving  the  two  incident  beams  to  add  vectorially.  In  Fig.  47, 
the  diagram  at  the  left  shows  the  momentum  of  the  collinear 
beams  adding  vectorially  to  the  sum  frequency  vector.  The 
same  Fig.  47  also  shows  the  vector  addition  of  the  momentum 
of  two  crossed  beams.  Notice  that  the  crossed  beams  will  not 
add  to  the  sum  frequency  vector  without  the  additional 
momentum  supplied  by  the  turbulent  eddy. 


FIG.  47.  Momentum  diagrams  for  collinear  and  crossed  beams 

In  the  collinear  case,  it  is  possible  to  add  the  incident 
wave  vectors  and  have  the  result  ks  =  (k0i  +  k02>  which  is  the 
wave  number  of  the  scattered  wave.  If  the  two  beams  cross  at 
an  angle  6,  the  magnitudes  of  the  incident  beams  add  to 

\X  +  ^  *  2k01k02  cos  6  ( A34) 


i 
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which  is  less  than  (ko  1  +  ko2)-  This  is  true  for  any  angle  other 
than  180'.  Thus,  the  momentum  of  the  crossed  incident  beams 
alone  is  insufficient  to  create  the  sum  frequency  component. 

The  additional  momentum  needed  to  radiate  the  sum 
frequency  component  comes  from  turbulent  momenta 
excitation. 


APPENDIX  B:  FOURIER  COEFFICIENTS 
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The  finite  Fourier  series  for  the  Kx  and  Ky  terms  are 
presented  on  the  next  pages  in  tabular  format.  Each  group  of 
coefficients  is  calculated  by  the  method  outlined  in  the  spectral 
theory  (eqs.  21  to  24). 

The  following  definitions  are  used  throughout  these 
tables: 

Kx(0)  =  ko+  F(0)  where  F(0)  =  1  -  ^  cos  0  -  7  sin  0  , 

Ky(0)  =  G(0)  where  G(0)  =  -  ^  sin  0  -  7  ^  cos  0  , 

and  7  =  ko-/ko+.  (Bl) 

The  coefficients  are  used  in  the  following  equations  of  average 
frequency,  variance,  third  moment,  and  fourth  moment  to 
determine  the  turbulent  flow  velocities. 


Average  frequency 

<  (0(1  >  -  0)04-  =  -  KjV ox  *  KyV oy 

=  -  ko+  F(0)Vox  -  ko+  G  (0) Voy  (B2) 


Second  Moment 

o2  =  <(g* i  -  <  oid  >)2>  =  ko+2[F2(0)<Vx2>  +  G2(0)<Vy2>  + 

2F(0)G(0KVxVy>]  (B3) 


Third  Moment 

<(o>d  -  <  cod  >)3>  =  -  ko+3[F3(0)<Vx3>  +  3F2( 0)G( 0)<Vx2Vy>  + 

3F(0)G2(0)<VxVy2>  +  G3(0)<Vy3>]  (B4) 

Fourth  Moment 

<(cod  -  <  tod  >y*>  =  -  ko+4[F4(0)<Vx4>  +  4F3(0)G(0KVx3Vy>  + 

6F2(0)G2(  0)<  V  X2V  y2>  +  4F(  0)G3(  0)<  Vx3  Vy> 
+  G4(0)<Vy4>]  (B5) 


In  the  following  Tables  (I  -  IV)  the  short  hand  notation  GG 
G2(8),  FG  =  F(0)G(0),  etc.  is  used. 
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1.  Correlator  -  Macintosh  program  to  compile  a  list  of 

velocity  correlations 

2.  Mean  velocity  correlations  -  Macintosh  program  that 

uses  the  least  squares  method  to  find  the  axial  and 
mean  velocities. 

3.  Rotator  -  Apple  He  program  that  automates  the 

rotational  experiment. 

4.  Analyzer  -  Apple  lie  program  to  calculate  spectral 

moments. 

5.  Mac  -  Pascal  program  used  to  calculate  spectral 

moments. 
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CORRELATOR  PROGRAM 
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'This  program  compiles  the  velocity  information 
’from  one  rotational  data  run  into  a  single  text  file. 

'Also  computes  the  rms,  skewness,  and  kurtosis 
'values  of  the  turbulent  parameters 

i 

DEFDBL  a-z 

DIM  variance(3),  mean(2),  cubic(4),  kurt(5)  'input  functions 
DIM  rms(3),  rmsMean(4),  S(4),  K(5)  'output  functions 

GO  SUB  inputTitle 

PRINT  "Will  you  want  to  output  to  disk  (y/n)? 

INPUT  "”;ans$ 

IF  NOT  LEFT$(ans$,l)  ="n"  THEN  GOSUB  outfile  ELSE  nofile  =  1 

GOSUB  readcorrelations 

GOSUB  outputRms 

GOSUB  outputRmsmean 

GOSUB  outputSkewness 

GOSUB  outputKurtosis 

INPUT"press  return  to  quit";a$ 

IF  nofile  =  0  THEN  CLOSE  #2 
END 

outfile: 

outf$  =  title$+":run_"+runNo$+":run_"+runNo$+". correlations" 
OPEN  outf$  FOR  OUTPUT  AS  #2 

nofile  =  0  ’boolean  variable  used  to  determine  whether  or  not 

to  save  to  disk 

RETURN 

inputTitle: 

'subprogram  to  input  data 
PRINT"What  is  tne  disk  title?"; 

INPUT  title$ 

PRINT  "  Run  #  =  ? 

INPUT  runNoS 


RETURN 


readcorrelations: 
subDir$(l)  =  "avg  Freq  Data" 
subDir$(2)  =  "sigma  Data" 
subDir$(3)  =  "skewness  Data" 
subDir$(4)  =  "kurtosis  Data" 

F0Rj%  =  2T0  5 

f$  s  "trident  2:"+subDir$(j%-l)+":run_"+runNo$+".cor" 
PRINT  "filename  =  ",f$ 

OPEN  f$  FOR  INPUT  AS  #1 
FOR  i%  =  1  TOj% 

INPUT  #1, label$  ’read  variable  descriptor 
INPUT  #1,  a$  'read  value 
IF  nofile  =  0  THEN 
PRINT  #2,label$ 

PRINT  #2,a$ 

END  IF 
a  =  VAL(a$) 

PRINT"  value  of  a  =\a 
SELECT  CASE  j% 

CASE  2 

mean(i%)  =  a 
CASE  3 

variance(i%)  =  a 
CASE  4 

cubic(i%)  =  a 
CASE  5 

kurt(i%)  =  a 
END  SELECT 
NEXT  i% 

CLOSE  #1 

IF  nofile  -  0  THEN  PRINT  #2,CHR$(13)  'Blank  line 
NEXTj% 

RETURN 


1  0  8 

outputRms: 

FOR  i%  =  1  TO  3 
SELECT  CASE  i% 

CASE  1 
type$=”xx" 

CASE  2 
type$  ="yy" 

CASE  3 
type$  =  "xy" 

END  SELECT 
IF  variance(i%)  <  0  THEN 
rms(i%)  =  SQR(-l*variance(i%)) 

ELSE 

rms(i%)  =  SQR(variance(i%)) 

END  IF 

PRINT  "RMS";type$;"  =",rms(i%) 

IF  nofile  =  0  THEN  'save 
PRINT  #2,  "RMS";type$ 

PRINT  #2,  rms(i%) 

END  IF 
NEXTi% 

PRINT 

IF  nofile  =  0  THEN  PRINT  #2,CHR$(13)  'Blank  line 
RETURN 

outputRmsmean: 

PRINT  "RMSxx/MeanX  =  \rms(l)/mean(l) 

IF  nofile  =  0  THEN  'save 
PRINT  #2,  "RMSxx/MeanX" 

PRINT  #2,  rms(l)/mean(l) 

END  IF 

FOR  i%  =  1  TO  3 
SELECT  CASE  i% 

CASE  1 
type$="xx” 

CASE  2 


type$  ="yy" 

CASE  3 
typeS  =  ”xy" 

END  SELECT 

PRINT  "RMS";type$;"/MeanY  =”,rms(i%)/mean(2) 
IF  nofile  =  0  THEN  save 
PRINT  #2,  "RMS";type$;7NifanY  " 

PRINT  #2,  rms(i%)/mean(2) 

END  IF 
NEXTi% 

IF  nofile  =  0  THEN  PRINT  #2,CHR$(13)  'Blank  line 

PRINT 

RETURN 


outputSkewness: 

FOR  i%  =  1  TO  4 
numerator  =  cubic(i%) 
denomX  =  (SQR(variance(l)))A(4-i%) 

-  denomY  =  (SQR(variance(2)))A(i%-l) 
SELECT  CASE  i% 

CASE  1 
type$="xxx" 

CASE  2 
typeS  =”xxy" 

CASE  3 

typeS  =  "xyy" 

CASE  4 

typeS  =  "yyy" 

END  SELECT 

skw  =  numerator/(denomX  *denom  Y) 
PRINT  "S";type$;"  =",skw 
IF  nofile  =  0  THEN  ’save 
PRINT  #2,  "S";type$ 

PRINT  #2,  skw 
END  IF 
NEXTi% 
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IF  nofile  =  0  THEN  PRINT  #2,CHR$(13)  ’Blank  line 

PRINT 

RETURN 

outputKurtosis: 

FOR  i%  =  1  TO  5 
numerator  =  kurt(i%) 
denomX  =  (SQR(variance(l)))A(5-i%) 
denomY  =  (SQR(variance(2)))A(i%-l) 

SELECT  CASE  i% 

CASE  1 

type$="xxxx" 

CASE  2 

type$  ="xxxy" 

CASE  3 

type$  =  "xxyy" 

CASE  4 

type$  =  "xyyy" 

CASE  5 

type$  =  "yyyy" 

END  SELECT 

krt  =  numerator/(denomX*denomY) 

PRINT  "K";type$;"  =",krt 
IF  nofile  =  0  THEN  'save 
PRINT  #2,  "K";  type$ 

PRINT  #2,  krt 
END  IF 
NEXTi% 

RETURN 


MEAN  VELOCITY  CORRELATIONS 


REM  Program  mean  velocity  correlations 
REM  This  program  computes  a  linear  square  fit  for  a,  b 
REM  which  are  the  time  averages  of  the  velocity  correlations 
xMean,  yMean 

'Note:  other  3  programs  to  calculate  velocity  correlations  using 
the  2nd,  3rd,  and  4th 

'moments  are  identical  to  this  one  except  that  the  functions  of 
Kx  and  Ky  are 

'slightly  more  complex,  the  sums  created  in  defineMatrix  are 
'different  but  follow  the  least  linear  squares  theory, 

’and  the  G-J  routine  has  to  reduce  larger  matrices.  But  since 
that  flexiblity 

'has  already  been  built  into  the  basic  structure  of  the  program, 

they  have  not  been 

'documented  with  as  much  detail. 

*************** 

'Define  variables  and  constants  for  equations 

theta  =  array  of  angles  at  which  Doppler  shifts  are  found 
DopplerFreq  =  array  of  time  averaged  Doppler  shifts  of 
scattered  spectra 

b  =  2D  array  for  linear  regression  matrix 
c  =  array  for  linear  regression  answers 
fminus  and  fsum  =  frequencies  (in  Hz)  of  the  difference  and 
sum  frequencies 

1482  =  speed  of  sound  (in  m/s)  in  water 
kfactor  =  constant  removed  from  b  array  to  speed  up 
calculations 
DEFDBL  B-c 

DIM  theta(lOO),  B(3,3),  c(3),DopplerFreq(100) 

LET  pi  =  3.1415927# 

LET  lastRow%  =  2 
LET  lastCol%  =  2 
LET  fminus  =  93665& 


LET  fsum  =  4107553& 

LET  factor  =  fminus/fsum 

LET  kfactor  =  fsum/  1482  '1482  m/s  =  speed  of  sound  in  water 

*********** ***** 

'Functions 

'  these  functions  create  the  theoretical  curves  that  the 

Doppler  shifts  follow 

DEF  FNConvertRadian#(x)  =  x  *  pi/  180 

DEF  FNKx#(rad)  =  -(1  -  SQR(2)/2*COS(rad)  -  factor  *  SIN(rad)) 
DEF  FNKy#(rad)  =  -(-SQR(2)/2*SIN(rad)  -  factor  * 
SQR(2)/2*COS(rad» 

***************** 

'Main  program 

Loads  the  data  into  memory 
sets  b  &  c  matrices  to  zero 

Defines  the  linear  regression  Matrix  as  groups  of  sums 

Prints  the  unsolved  matrices 

Uses  Gauss-Jordan  row  reduction  to  solve  matrix 

Prints  the  solve  matrix 

Prints  the  velocity  solutions 

Creates  a  data  file  to  compare  regression  values  to  actual 

data 

GOSUB  loadSigma 
FOR  i  =  1  TO  lastRow% 

FOR  J  =  1  TO  lastCol% 

LET  B(i,  J)  =  0 
NEXT  J 
LET  c(i)  =  0 
NEXT  i 

GOSUB  defineMatrix 

GOSUB  printMatrix 

CALL  gaussReduc  (lastRow%,  B(),  c()) 

GOSUB  printMatrix 
GOSUB  outputAnswer 
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GOSUB  errorAnalysis 
END 

<****************** 

’sub  program  to  input  data  from  files 

’  program  assumes  that  zenith  analyzed  the  data  and  created 
'  the  file  under  subdirectory  ,,RUN_#"  and  file  name 
"RUN_#.avg” 

'  where  the  #  is  a  number  1-14 
'  Variables  affected: 

'  theta(),  DopplerFreqO  have  proper  values  for  the  run# 

'  N  =  number  of  angles  for  that  run 
loadSigma: 

PRINT’What  is  the  disk  title?"; 

INPUT  title$ 

PRINT  "  Run  #  -  ? 

INPUT  runNo$ 

f$=  title$+":RUN_,’-HTjnNo$+":RUN_"+runNo$+".avg" 

PRINT  "filename  =  ";f$ 

OPEN  f$  FOR  INPUT  AS  #1 
k=0 

'Note:  each  line  of  data  looks  like 
'  angle  <tab>  doppler  shift  <CR> 

’  so  the  parser  routine  below  simply  finds  the  <tab>  =  chr$(9) 

'  and  throws  the  left  half  in  a$,  and  the  right  half  in  b$ 

WHILE  NOT  FOF(l) 

INPUT  #l,a$ 
i%=  0 
done  =  0 

WHILE  (i%<  LEN  (a$))  AND  (done  =  0) 
i%  =  i%+l 

IF  MID$(a$,  i%,l)  =  CHR$(9)  THEN  done  =  1 
WEND 

a=  VAL(LEFT$(a$,i%)) 

B$=  RIGHT$(a$,  LEN(a$)-i%) 

B  =  VAL(B$) 
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IF  (LEN(a$)  o  0)  THEN 
theta(k)  =  a 

DopplerFreq(k)  =  B  -  fsum 
k=k  +  1 
END  IF 
WEND 
CLOSE  #1 
N  =  k  -  1 
RETURN 

'**************** 

'sub  program  to  output  the  mean  velocities  after  the  gauss- 
'Jordan  reduction  of  the  matrix  is  complete 
'Since  G-J  only  explicitly  solves  for  the  last  variable,  the 
'program  uses  the  lastRow  to  solve  for  the  variable  above  it. 

Then  it 

'cascades  upward  using  solutions  that  have  already  been  found 
to 

'get  the  next 

'Also  before  printing  a  kfactor  is  removed  from  the  solution, 
this  factor  was 

'removed  from  the  original  matrix  B  to  avoid  having  the 

computer  create  overflow 

'errors. 

'Variable  affected 

c()  =  velocity  correlations  upon  exit 

outputAnswer: 
coeff$(l)  ="xMean" 
coeff$(2)  ="yMean" 

PRINT"From  gaussian  reduction  the  solution  to  matrix  equation 
is" 

c(lastRow%)  =  c(lastRow%)/  B(lastCol%,  lastRow%) 

FOR  row  =  (lastRow%  -1)  TO  1  STEP  -1 
sum  =  0 

FOR  col  =  row+1  TO  lastCol% 
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sum  =  sum  +  B(row,  col)*c(col) 

NEXT  col 

c(row)  =  c(row)  -  sum 
NEXT  row 

FOR  row  =  1  TO  lastRow% 

PRINT"  _ " 

c(row)  =  c(row)  /  (kfactor) 

PRINT" coefficient  ";coeff$(row);"  =  ”;c(row) 

NEXT  row 
PRINT 

'This  last  section  allows  the  user  to  save  the  calculated 
’velocity  information  to  a  text  file  along  with  a  label  -Coeff$() 

'for  each  velocity 
'file  format  is: 

'name  of  velocity  correlation  <CR> 

'value  of  velocity  correlation  <CR> 

'  ...repeat  for  all  velocities 

PRINT  "Do  you  want  to  save  to  disk?" 

INPUT"  ";ans$ 

IF  NOT  LEFT$(ans$,  1)  ="n"  THEN 
PRINT  "Default  saves  to  DISK  =  trident  2" 

PRINT  "  SubDirectory  =  Avg  Freq  Data" 

PRINT  "  filename  =  run_";runNo$;".cor" 

PRINT 

PRINT  "Do  you  want  to  use  the  default  save? 

INPUT  "”;default$ 

ELSE 

GOTO  endOutputAnswer 
END  IF 

IF  LEFT$(default$,l )  ="n"  THEN 
f$  =  FTLES$(0) 

ELSE 

f$  =  "trident  2:Avg  Freq  Data:run_"+runNo$+".cor" 

END  IF 

OPEN  f$  FOR  OUTPUT  AS  #1 
FOR  i  =  1  TO  lastRow% 
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PRINT  #1,  coeff$(i) 

PRINT  #1,  c(i) 

NEXT  i 
CLOSE  #1 

endOutput  Answer: 

RETURN 

defineMatrix: 

***************** 

'sub  program  to  initialize  the  b  &  c  matrices 

'Each  element  of  b  &  c  are  defined  by  a  sum  of  function  values 

from 

'the  linear  least  squared  fit  theory. 

'This  routine  sets  up  the  appropriate  functions  and  sums  them 
over 

'all  N  data  points  for  the  run. 

'Variables  affected 

b(),  and  c()  which  are  set  to  initial  values 
FOR  loop  =  1  TO  N 

angle  =  FNConvertRadian#(theta(loop)) 
c(l)  =  c(l)  +  FNKx#(  ang  le)  *  DopplerFreq(  loop) 
c(2)  =  c(2)  +  FNKy#(angle)*DopplerFreq(Ioop) 

B(l,l)  =  B(l,l)  +  FNKx#( angle)*  FNKx#(angle) 

B(l,2)  =  B(l,2)  +  FNKx#(angle)*  FNKy#(angle) 

B(2,2)  =  B(2,2)  +  FNKy#(angle)*  FNKy#(angle) 

NEXT  loop 

'Note:  b  is  a  symmetric  matrix,  thus  values  repeat  in  the  matrix. 

'To  avoid  having  the  computer  calculate  the  sum  again  these 
values  are  simply 

'transcribed  after  the  sum  is  calculated  once. 

B(2,l)  =  B(l,2) 

RETURN 


printMatrix: 

****************** 
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'sub  program  to  output  the  current  values  of  the  matrix  B 
The  matrix  is  output  for  debugging  purposes  and  to  insure 
'the  validity  of  the  G-J  reduction 
"Variables  affected 
'  none 

FOR  row  =  1  TO  lastRow% 

FOR  col  =  1  TO  lastCol% 

PRINT  B(row,  col);" 

NEXT  col 
PRINT  c(row) 

NEXT  row 

PRINT 

RETURN 

SUB  switch(x%,  y%,  B(2),  c(l),  last%)  STATIC 
<****************** 

'sub  program  to  interchange  rows  x,  y 

'if  a  zero  were  to  occur  in  the  matrix  during  the  G-J  reduction, 

'it  could  cause  a  division  by  zero  error.  This  routine  switches 
out  such 

'a  row  for  a  non-zero  one. 

'Variables  affected 

b(),  c()  which  have  2  rows  interchanged 
FOR  col  =  1  TO  last% 
q#  =  B(x%,  col) 

B(x%,  col)  =  B(y%,  col) 

B(y%,  col)  =  q# 

NEXT  col 
q#=  c(x%) 
c(x%)  =  c(y%) 
c(y%)  =  q# 

END  SUB 

errorAnalysis: 

******************* 

'This  program  calculates  the  best  fitted  values  from  the 
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'computed  velocity  correlations.  It  then  saves  the  best  fit 
values  and 

‘the  actual  data  to  a  text  file  so  the  user  can  examine  the 
accuracy  of  the  fit. 

'Data  file  has  the  form 
’*  {Data  headings}  <CR> 

’angle  at  which  data  taken  <TAB>  Regression  value  <TAB>  Actual 
data  <CR> 

’....repeated  for  all  N  data  points 
fout$  =title$+":RuN_’'+ranNo$+":errorA" 

OPEN  fout$  FOR  OUTPUT  AS  #3 
'Make  data  headings  for  Cricket  graph 
'Data  headings  start  with  a  * 

PRINT  #3, 

outst$  =”theta  star"+CHR$(9)+"Regression  Plot"+CHR$(9)+"  Actual 
Data” 

PRINT  #3,  outst$ 

'Since  kfactor  has  been  removed  from  the  fitted  constants,  it 
must  be  factored 
'back  in. 

FOR  loop%  =  1  TO  N 

angle  =  FNConvertRadian(theta(loop%)) 
f  =  +c(l)  *FNKx#(angle)  +  c(2)  *  FNKy#(angle) 
f  =  f*kfactor 
q  =  DopplerFreq(loop%) 
outst$  =  STR$(theta(loop%))+CHR$(9)  + 

STR$(f)+CHR*(9)+STR$(q) 

PRINT  #3,outst$ 

NEXT  loop% 

PRINT 
CLOSE  #3 
RETURN 


SUB  FindZero  (topRow%,  B(2),  c(l),  errors$,  lastRow%)  STATIC 

***************** 


'sub  program  that  checks  for  zeroes  in  the  lead  row 
‘calls  switch  if  necessary  to  avoid  division  by  zero  errors. 
’Variables  affected 

errorsS  is  a  boolean  string 

=  'false'  if  no  zeros  found  or  the  row  was  replaceable 
=  'true'  if  could  not  replace  row  which  occurs  if  the  matrix 
is  nonsingular  (i.e.  no  unique  solution) 
b(),  and  c()  could  have  2  rows  interchanged 

col%  =  topRow% 
errors$=  "false" 

WHILE  (B(topRow%,  col%)  =  0)  AND  (  errors$  =  false") 
row%  =  topRow%  +  1 
check$  ="working" 

'the  lead  row  had  a  zero,  so  now  we  must  look  through  the 
remaining  rows  to  find 

'an  non-zero  lead  row  to  replace  it  with. 

WHILE  checks  =  "working" 

IF  (B(row%,  col%)  =  0)  THEN  row%  =  row%  +  1  ELSE  checks  = 
"done" 

IF  row%  >  lastRow%  THEN  checks  =  "done” 

WEND 

IF  NOT  row%  >  lastRow%  THEN 
switch  row%,  topRow%,  B(),  c(),  lastRow% 

ELSE 

errors$=  "true" 

END  IF 
WEND 
END  SUB 


SUB  gaussReduc  (N%,B(2),  c(l))  STATIC 
***************** 

’sub  program  to  perform  G-J  reduction  upon  a  NxN  matrix 
'and  its  solution  vector 

'Note:  in  this  case,  B  is  2  dimensional  but  the  algorhythm  works 


'for  any  square  matrix 

’Procedure:  make  first  column  a  1  and  zeros 

then  move  to  second  column,  second  row  and  do  the 

same, 

cascading  downwards 
'Variables  affected: 

'  b()  is  half  diagonal 

c()  is  manipulated  to  allow  b()  to  be  half  diagonal 
nonsing$  is  used  to  talk  between  findZero  and  gaussReduc 
to 

detect  nonsingular  matrices. 

FOR  row  =  1  TO  (N%-1) 

leftCol%  =  row  ‘use  leftRow  simply  for  readability 
findZero  leftCol%,  B(),  c(),  nonsing$,  N% 

IF  nonsing$="true"  THEN  breakOut 
invers#  =  1/  B(row  ,  leftCol%) 

B(row,  leftCol%)  =  1 
FOR  col  =  leftCol%+l  TO  N% 

B(row  ,  col)  ss  B(row  ,  col)  *invers# 

NEXT  col 

c(row)  =  c(row)  *  invers# 

FOR  nextRow  =  row+1  TO  N% 

multiplier#^  -l*B(nextRow,  leftCol%) 

FOR  col  =  leftCol%  TO  N% 

B(nextRow,  col)  =  B(nextRow,  col)  +  B(row,  col)*multiplier# 
NEXT  col 

c(nextRow)  =  c(nextRow)  +  c(row)  *  multiplier# 

NEXT  nextRow 
NEXT  row 
breakOut: 

IF  nonsing$  =  "true"  THEN  PRINT  "Oh  NO!!  nonsingular  matrix 
cannot  row  reduce!!!!" 

END  SUB 


ROTATOR  PROGRAM 
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1  REM  **•*•**. ...... 

£  REM  ROTATOR.  £.  MAC.  -oR  IQhC'G 
4  REM  MODIFIED  03APR3E  •  • 

6  REM  LINES  1  —  100  RESET  AND  RECONFI GUR:":  A*.,.  SIG--L__-  l.  -  -I'-L-'"'  ► 

PPLE  J ie 

IE*  REM  ALPHA*  3L'Nf AinS  SPECIAL  KEYSTROKES  FOR 

1£  REM  MEMORY  LOCATIONS  -16£63  TO  -16E36  ARE  DSUiUl'.t-  .  . 

IATOR  SWI TCHES  0-3.  PEEnmc  THE  EVEN  VAlUES  TU-v.'i:'  :  ml;’!  OF: 

.  ODD  TURNS  THEM  UN. 

1 H  REm  LINES  46  —  £4  SET  UP  SCAN  LOOP'S  AND  TIMING  FOR  STEPPER  MOTOR. 


26 

4£ 

44 

46 

■46 


50  gi 


DIM  ALPHA*  l  £f  )  :  GOSUB  6900:  REM  8300  -  INITIALIZE  ALPHA* 


l6£90>  :A  =  t-'EEK  (  -  16£3£)  :( 


REM  GET  INITIAL  DATA 


:'EcK  i  -  j6E34j  :A  =  PEEI 


val 


54 

56 

56 

6« 

6£ 

64 

66 

66 

7d 

7£ 

74 

76 

76 

60 

35 

36 

37 


38 


A  =  PEEK  (  • 

<  --  16£36; 

GOSUB  6500: 

HOME 

A»  *  "number  of  7L5  sweeps  per  angle  “:LB*  =  GOSUB  7000:NS 

(MM* ) 

GOSUB  6£00:  REM  GET  INI T I AL, F INAL,  AND  ANGULAR  INCREMENT. 

NA  «  INT  <<AF  -  AI)  /  AD):NA  =  NA  +  1 
PRINT  "number  of  angles*  ";NA 

—  O  /  4)  *  (4  /  6)  *  (600)  *  (£00):  PRINT  ST ;  "steDS=  1  deg  or*  h— be 

am  apparatus  ":ST  *=  ST  *  AD 

PRINT  "Number  of  steps  to  increment  ” ; AD; “  degrees  « 

INT  (ST):  PRINT  "integer  value  for  number  of  steps* 


AD  /  ST:  PRINT  "correspond i ng  to  actual  increment  angle (deg) 
INT  ( (MS  -  Ml  *  65536)  /  £56>:M3  *  MS  -  M 


MS 

AC  *  MS 
=  "  ;  AC 

Ml  -  INT  (MS  /  65536) :M£  - 
1  *  65536  *  M£  *  £56 
GOSUB  6300:  REM  GET  DATA 
GOSUB  15000:  REM  DECIDE 
GOSUB  6000:  REM  INPUT  PULSE 
GOSUB  1 0£00 :  REM  PRINT 


;ST 
'  ;MS 


D*  = 
Z*  « 
I*  «= 
J*  * 
REM 
REM 
REM 


CHR*  (4) 
CHR*  (£6) 
CHR*  (3£ 
CHR*  (3£ 


TITLE 

WHETHER  OR  NOT  TO  SAVE  TRIGGER 
LENGTH 

ACTUAL  MICROPULSE  LENGTH 


RAMPS 


I  ) 
3£ 


I ) 


LINES  100-540  REPEATED 


FOR  EVERY  ANGLE  SCANNED 

DELAY  LOOPS  ARE  TO  ALLOW  IQ400  TIME  TO  EXECUTE  COMMANDS  BEHOOF 
APPLE  SENDS  ANOTHER  COMMAND.  WITHOUT  THE  DELAYS.  ThF  lk'jGL 

AM  WILL  CRASH' 


REM  1300  -  TURNS  ON  IEEE 
"WT"  -  COMMAND  TO  IEEE 
APPLE 


1600  -  TURNS  OFF  I  Eft 
‘ RD"  -  IEEE  OJTPUT  TO 


33  REM 

100  FOR  Q  =  1  TO  NA  STEP  1 
105  AN  *  A I  +  ( O  -  1)  *  AD 

107  GOSUB  3000:  REM  CREATE  IQ400  MACROS 

108  REM  LINES  110-175  RE-ZERO  ID400  ARRAYS  BETWEEN  ANGLES 
110  GOSUB  1300 

140  PRINT  "WT" ; I*; 2*;Al* 

150  PRINT  "WT"; I*;Z*;B1* 

160  PRINT  "WT";1*;Z*;  A£i> 

170  PRINT  "WT" ; 1*;Z*;B£* 

175  GOSUB  1600 
178  REM 


12- 


r£m  lines  160  -  330 

FOR  K  =  1  TO  NS  ST  tP  1 

REM  LINES  IBS  -  330 

D  IQ400 

GOSuB  1300 

PRINT  '  Wf  ;  I*:Z*:"R" 

GOSuB  l60i- 
GOSUB  4000 

REM  time  ue  *  a v  =  uUsec 
FOR  J  =  i  TO  39000  S'  EP  1  : 
GOSUB  1300 

PRINT  :  RE’! 

GOSUB  1  600 

REM  aeiay=  3  sec 

FOR  J  =  1  TO  5000  STEP'  1: 

REM  LINES  340  -  360 

GOSUB  1300 

PRINT  "WT" ; IS:Z*:A3* 

GOSUB  1600 
REM  de 1 ay=3sec 


REPEATED  FOR  EVEKi  7l5  SCAN 
COORDINATE  TRIGGERS  BETWEEN  THE  7L5  AN 


HE  HOLD  BUTTON 


;  1 LL  TRIGGER 


NEXT  j 

SIGNAl. 


AVERAut 


-i c.  E  IC'^00  ARRAv  i> 


34B 

FOR  J 

=  1  TO  8000  STEP  1 : 

NEXT  J 

ESC 

IF  NOADD  =  1  THEN  390:  h£M 

IF  v 

i  oGER 

RAMPS  THEN  SKIP 

353 

GuSLi £< 

1300 

35** 

PR  I  N*i 

" wT " : 1  S : Z*  : B3* 

356 

GOSU& 

1600 

358 

REM  deiay=3sec 

360 

FOR  J 

=  1  TO  6000  STEP  1  : 

nE XT  J 

390 

PRINT 

"10400  LOOP=  f K 

300 

PRINT 

"APPLE  lie  LOOP'=  "  ;  0 

310 

PRINT 

"SCATTERING  ANGLE = ' : 

330 

next  i 

K 

333 

rt£ M  l.  I NE5  ^35  —  A&v? 

S^vEi 

335 

GOSUB 

1300 

34C 

PRINT 

"WT" ; IS;Z*:A4* 

345 

GOSUB 

1  600 

35<? 

REM  delay=lsec 

360 

FOR  J 

s  i  TO  3600  STEP'  1  : 

NEXT  J 

365 

GOSUB 

1 300 

370 

PRINT 

"WT" ; 1 *; ZS ;B4» 

375 

GOSUB 

1600 

380 

REM  delay=lsec 

390 

FOR  J 

=  1  TO  3600  STEP  1  : 

next  j 

395 

GOSUB 

1  300 

400 

PRINT 

"WT" ; I*;Z*;ClS 

4C3 

GOSUB 

1600 

405 

FOR  J 

=  1  TO  3000:  NEXT  J 

430 

FOR  J 

=  1  TO  3000  STEP  1  : 

NEXT  j 

435 

GOSUB 

1300 

440 

PRINT 

"WT" ; I* : Z* ;«5» 

445 

GOSUB 

1600 

450 

REM 

de 1 ay=4sec 

460 

FOR  J 

=  1  TO  15000  STEP  i  : 

ne  ».  r  j 

461 

GOSUB 

5000:  REM  HARD  RESET  SENT 

CONTROL  AFTER  DISK  USt 

480 

REM 

delay=4sec 

490 

FOR  J 

=  1  TO  15000  STEP  1 : 

:  NExT  J 

500 

REM 

a  space  to  later  or  ir.t  out  a 

510 

REM 

a  space  to  put  i r  an 

60s ec  at 

530 

IF  0 

<  ( NA  >  THEN  GOSUB  6000:  REn' 

THE  TR 


SAvES  AVERAGED  SPtCTRA  TO 


REGAIN 


6000  -  ADVANCES  STEPPER 


MOTOR 


125 


ANGLE 


t»40 

900 

910 

1290 

1295 

1297 

1300 

1310 

1320 

1330 

1340 

1490 

1495 

1497 

1500 

1510 

1520 

1590 

1595 

1597 

1600 

1620 

1630 

1640 


TO  NEXT 
NEXT  0 
GOSUB  5000 
END 

REM  ******************** 

REM  1300  SUBROUTINE  TURNS  ON  IEEE  CARL 
REM  ******************** 

PRIN'  D*:"PR#7" 

PRINT  D*;"IN#7" 

PRINT  ••SCI" 

PRINT  “RA" 

RETURN 

REM  ******************** 

REM  1500  SUBROUTINE  RESETS  APPLE  PERIPHERALS 

REM  *****••**»••**••*•»* 

PRINT  D«;“IN«0" 

PRINT  D*;"PR#0" 

RETURN 

REM  ****»*•**••»*»«••••• 

REM  1600  SUBROUTINE  TURNS  OFF  IEEE  CARD 

REM  •*•*•«••**••*•••*•*• 


PRINT  -UT" 

PRINT  D* ; “ IN#0" 
PRINT  D*;"PR#0" 
RETURN 


2990 

REM 

2995 

REM 

3000  SUBROUTINE 

SETS  UP  1 0400  MA, 

2997 

REM 

3000 

REM 

Strings  cf  ASCII  code  to 

3010 

REM 

remot  e  cont  rc.  1 

the  10400 

3020 

Al*  - 

"  "HAGAI  •' 

3030 

Bl*  * 

" -KCGCI" 

3040 

A2*  = 

"  A60-AIAI" 

3050 

B2*  = 

AB0=CICI ” 

3060 

A3*  *> 

'  AiiAGAIAI" 

3070 

B3*  * 

"  AwCGCICI" 

3080 

A4*  = 

"  ACPI"  ♦  STR* 

(NS)  ♦  "-AI" 

3090 

B4*  = 

"ACC I"  ♦  STR* 

(NS)  ♦  "*CI" 

3100 

G2*  = 

STR*  (AN) 

3102 

Gl*  * 

"  K" 

3104 

G3*  = 

•=M00" 

3106 

Cl*  * 

Gl  *  +  Gel  $  -*•  G3S 

3110 

A5*  * 

" _AAC"  +  NAM*  + 

"-M00AAI " 

31 30 

C3*  = 

"  AAC"  ♦  NAM*  + 

"“M00EC01 = " 

3140 

C4»  = 

"  DAE" 

3150 

RETURN 

3990 

REM 

►  ***♦■* 

3995 

REM 

4000  SUBROUTINE 

TURNS  ANNUNCIATOR 

R  7L5’ 

S  SWEEP 

4000 

REM 

4010 

A  «  PEEK  (  -  16293) 

FOR  J 

*  1  TO  50:  NEXT 

J 

4030 

A  =  PEEK  <  -  16294) 

RETURN 

4990 

REM 

***** 

4995 

REM 

5000  SUBROUTINE 

TURNS  ANNUNCIATOR 

ESET  T 

G 

5000 

REM 

1  On  X  OFF  TO 


X  OFF  to  send 


5010 

5020 

5030 

5035 

5040 


A  = 
FOR 
A  * 
FOR 


PEEK 
J  =  1 
PEEK  ( 


<  -  1 629 1  > 

TO  50:  NEXT  J 
-  16292) 


RETURN 


J  =  1  TO  4000:  NEXT  J 


TRIGGE 


HARD  R 


EXECUTES  A  MACHINE  LANGUAGE 


5995  REM  6000  SUBROUTINE 

E  TO  ACTIVATE  STEALER  MOTOR 
6000  REM  SUBROUTINE  IS  CALLED: 

6010  REM  MSK. FEB19S1.  STEPPER 
6060  REM  PUlSER  PROGRAM 
6030  REM  FOR  USE  W/  MOTOR  100  ON  1500  OFF 
6040  REM  Si_0-SvNC  STEF'  MOTOR  WORKED 

6050  REM  WELL  t*  50M 1 CROSEC " ON"  AND  1 500M I CROSEC " OFF " 

6060  REM  INPUT  f  MICROSECONDS  ON 
6070  REM  INPUT  F  MICROSECONDS  uFF 

6060  rem  number  of  bursts= 

6090  REM  S  i  EPS  ARE  _  I  Ml  TED  TO  656-*- 1 656  '  •  i  656  > 

6100  REM  DAT£:FEB  6.  T'9'9i 

6il0  REM  USE  Annul  : r*  OR  PORT  *0  FOR  PUlSES 
6130  REM  THIS  IS  A  MODIFICATION  Or  TriE 
6l30  REMi  PULSE. BURS. . SEQUENCE  or'ogram 

6140  REM  PROGRAMMERS  =  MURRAY  kORMhN  »  ChRISTInA  BURCH 
6145  REM  ♦»*»»***«►*■■*■**♦•#■♦**♦♦* 

6150  GOSUB  10000:  RE  r  GET  J  AND  K  IlEnGTm  OF"  P'ulSESi 

6600  R  *  m3  : 5  =  m- 

6310  FOR  I  =  864  TO  69. 

6330  READ  A 
6330  POKE  I, A 
6340  NEXT  1 
6345  RESTORE 

6350  iF  R  —  0.  fnEX  ju  u  6*rl0 

6360  POKE  665. 1 

6370  POKE  667,  R 

6360  POKE  673.3 

6390  POKE  661, h 

6*t00  CALL  6fcP 

6K10  PR  IN  GE  •  -\6aDv 

6430  PhINi  ' M i Cm-jSElw  .t’s  On  =  :K 

6430  PRINT  MICROSECOND--  OFF  =  " :  J 

6n«0  PRINT  RcAOt 

6450  PRINT  GO.  .  . 

6460  IF  S  —  0  'iLEN  o_'  -  65**0 

6*. 70  POKE  667,  655 

6**60  POKE  865.  S 

6490  POKE  667, 655 

6500  POKE  673,6 

6510  POKE  661, H 

6530  CALL  664 

6540  IF  Ml  =  0  TrtEN  bo60 

6550  POKE  667, c55 

6560  POKE  665, 355 

6570  POKE  667 , 65- 

6580  POKE  873, G 

6590  POKE  BBl.R 

6600  CALL  864 

6610  Ml  *  Ml  —  i  :  GuT'j  65*** 

6660  FOR  W  =  1  TO  5:  i’Fu-  -  1 9ts :  NEX  .  w 
6630  P'RINT  Di  :  REM  bet  OF  Cal-L  **.  96 

6640  RETURN 

6990  REM 

6995  REM  7000  SUBROUTINE  OUTPUTSSTR I  NGS  10  -'  *<  1  N't  r .  -  ESSAR  r  BECAuS 

E  PRINTER  CARD  TURNED  UN  S  OFF  CUT  E  FREQUENT  . 

7000  REM  ********************** 

7010  PRINT  "  Enter-  :  A*  ;  ■  (";  LB* :  INPUT  MM1 


127 


7020  NN*  =  .  tf  14  iA»,  1):  REM  ALTERNATE  EN7  H  r  -'DInT  FOR  OTHER  RDLiT  i  isi _ 

7030  NN*  =  CHR*  <  ASC  INNS)  -  32) 

7040  A*  =  NN*  +  RIGHT*  (A*,  LEN  (A*)  -  1) 

7050  REM  LINES  7020-7040  CAPITALIZE  THE  FIRST  LETTER  OF  A* 

7060  PRINT  CHR*  (4>:"Pfi#l" 

7070  PRINT  '  -;A*s“  =  "  ;  NM* ;  "  "  ;LB* 

708i ZI  PRINT  CHR*  <4>;"PR#0" 

7090  RETURN 

7095  REM  *****♦«♦**♦♦♦*♦♦**♦** 

7097  REM  7100  SUBROUTINE  OUTPUTS  SIMPLE  STRINGS  TO  PRINTER 

7100  REM  ********************* 

7110  PRINT  CHR*  (4 )  ;  "PR# 1 11 
7120  PRINT  A* 

7130  PRINT  CHR*  (4);"PR#0" 

7140  RETURN 

7190  REM  ********************* 

7195  REM  7 £00  SUBROUTINE  OUTPUT  BLANK  LINE  TO  PRINTER 
7E00  REM  CHR* ( lO)  IS  A  LINE  FEED  CHARACTER 

7E05  REM  »*•••»**»••****••*•*• 

72  1 0  A*  -=  CHR*  (10) 

7££0  GOSUB  7100 
7£30  RETURN 

7£90  REM  ********************** 

7c' 95  REM  7300  SUBROUTINE  SETS  UP  OUTPUT  FOR  PRINTER  IF  TWO  OUTPUT 

S  ORE  NEEDED. 

7300  REM  ********************** 

7310  PRINT  *'  11  s  A 1  * ;  "  (  "  5 LBS 5  :  INPUT  Nit 

7320  PRINT  ••  "  5  A2* ;  "  (  "  ;LB*  INPUT  N£* 

7330  A*  -  ••  "  +  Alt: MM*  *  Nl*:  GOSUB  70SOs  REM  USING  ALTERNATE  ENTRY 

POINT 

7340  A*  «  "  "  +  A£*:MM*  *  N£»:  GOSUB  7060:  REM  USING  ALTERNATE  ENTRY 

POINT 

7350  RETURN 

7390  REM  ********************** 

7395  REM  7400  SUBROUTINE  SETS  UP  OUTPUT  FOR  CHOICES 

7400  REM  •***•*•«*•••*•••**•**• 

7*10  PRINT  "  Enter  ":A*;"  (  "  ; CH* ;  "  )  *'  s  :  INPUT  MM* 

7420  NN*  =  LEFT*  (A*, 1) 

7430  NN*  =  CHR*  (  ASC  <NN*>  -  3c) 

7440  A*  =  NN*  +  RIGHT*  <A*,  lEN  (A*>  -  1) 

7450  REM  LINES  7020-7040  CAPITALIZE  THE  FIRST  LETTER  OF  A* 

7460  PRINT  CHR*  <4>;"PR#1" 

7470  PRINT  M  " ; A* ; "  «  " ; MM* 

7460  PRINT  CHR*  (4>;"PR#0" 

7*90  RETURN 

6000  REM  ********************* 

6005  REM  8000  SUBROUTINE  INPUTS  PUcSE  SPEEDS  FOR  ROTATION  STEPPE 

R  MOTOR. 

6010  REM  ********************* 

8020  PRINT  :  PRINT  "default  sett  inns":  PAINT  for  the  Dulse  duration  a 

rid  pulse  repetition  period  of  the  electronic  pulse  to  the-  e 

tepper  motor.":  PRINT  “  200  microseconds  on.":  PRINT  "  6000  unerase 

conds  off 


INPUT  Is  this  acceptable  (y/r.) 


Not  e  : 


this 


wili  be  tne  last  prompt  before  the  experiment  commences.  ‘';D"~S 

S040  IF  LEFTS  (DFTS,  1)  *  "N"  OR  LEFT*  (DFTS,  1  )  =  "r,"  THEN  8060 

6050  0  =  800  :P  =  6*00:  GOTO  80S* 

8060  INPUT  "How  long  should  tne  Dulse  be  on  (m i crosec ) ?  Note  800  should 
be  the  minimum.  ";0 

8070  INPUT  "How  long  should  the  Dulse  be  off  'vote:  6000  should  be  mir, 
i  mum.  “ :  P 
8880  RETURN 

8! 90  REM  ***********■**■*■■»•*•■*■***•* 

8195  REM  6800  SUBROUTINE  GETS  3  IMPORTANT  PARAMETERS  AI  - 

INITIAL  ANGLE  OF  SCAN  AF  -  FINAL  ANGLE  OF  SCAN  AD  -  A 

NGULAR  INCREMENT  BETWEEN  SCANS 

8800  REM  ********************** 

8810  A*  =  "initial  angle":LB*  —  “deg'':  GQSUB  7000:AI  =  UAL  iMMS) 

8880  AS-  =  "final  angle":LB*  =  "oeg":  GOSUB  7000:AF  =  VAL  CMMS) 

8830  ?:■>•  =  "increment  angle":i_B*  =  "deg":  GOSUB  7000:AB  =  VAl  tMMS) 

8840  RETURN 

8895  REM  8300  SUBROUTINE  INPUTS  A  FILENAME  TO  SAVE  THE  DATA  U 

NDER. 

6897  LINESS380  -  8390C  ON  VERTTHEFI  LEN  AME  TO  ASTRI NGTHE IO400CANUSE. 

8300  REM  ********************** 

6305  NAMS  =  "" 

8310  PRINT  :  INP'uT  "Enter  a  5  character  code  to  save  the  data  under.  . 

.  ’T  ITLES 

8315  7-  LEN  (TITLES)  >  5  THEN  PRINT  “ERROR1 1  The  maximum  length  of  a 

data  label  is  5  characters.  GOTO  8305 

8317  IF  LEN  (TITLES)  <  1  THEN  PRINT  "ERROR11  You  must  tvoe  something' 

1  "  :  GOTO  8305 

8380  FOR  .  =  1  TO  LEN  (TITLES'. 

8330  CHARS  =  MIDS  (TITLES, J, 1 ) :  REM  GET  A  CHARACTER  OFF  OF  STRING 
8340  IF  (  BSC  (CHARS)  >  96)  AND  (  ASC  (CHARS)  <  183)  THEN  CHARS  =  Ct-tRS 
l  ASC  (CHARS)  -  38):  REM  MAKES  LOWER  CASE  INTO  UPPER  CASE 
6350  IF  (  ASC  (CHARS)  >  64 >  AND  (  ASC  (CHARS)  (  91)  THEN  NCHARS  =  ALPHAS 

i  ASC  (CHARS)  -  64):  GOTO  6380:  REM  IF  LETTER  THEN  DECODE 
6360  IF  (  ASC  (CHARS)  >  47)  AND  (  ASC  (CHARS)  (  56)  THEN  NCHARS  =  CHARS: 

GOTO  8380:  REM  A  NUMBER  HAS  BEEN  CHOOSEN 
6370  PRINT  :  PRIN’"  "ERROR1  1  You  can  only  enter'  a  1  Dhanumer  i  ccnar  act  er  s.  " 

:  GOTO  6305 

6380  NAMS  =  NAMS  ♦  NCHARS 
6390  NEX"  J 

8398  AS  =  "Data  file  name... ":Mms  =  TITLES:lBS  =  GOSUB  7060 

6395  RETURN 

6490  R£n  ********************** 

54  95  REM  8500  SUBROUTINE  GETS  7 1 ME  AnD  DATE  OF  SCAN.  ThEN  ASKS  I 

F  USER  WANTS  TO  INPUT  ALL  PARAMETERS  OF  THE  SCAN 
6500  RE*"  ********************** 

8510  INPUT  "Enter  date  ( Ju'N-08-  1 998 )  ";DTS 
5580  INPUT  "Enter  time  (HHMM)  ";TMS 
6585  PRINT  CHRS  (4>:"PR*i“ 

5530  PRINT  “Date  :  ’ :  D'  s  :  PRINT  l me :  "  ;  TMS 

6535  PRINT  CARS  (  4  )  ;  11  PR«0" 

65*0  INPUT  "Do  you  wish  to  enter  the  settings  for  the  7L5  spectrum  ana 
lyzer  and  the  i_eCroyn’  ;  ANS :  IF  LEFTS  (ANS,  1  )  =  ”Y"  THEN  GOSUB 

8600 

8550  RETURN 

6590  ROM  *********** ******** ♦  * 

af,  95  REM  8600  SUBROUTINE  GETS  SET  ; NGb  OF  MONITORING  EDO  I  PM 

ENT  OF  THE  APPARATUS 
6600  REM,  ********************* 

8610  GOSUB  7800: AS  *  “7l5  Spectrum  Analyzer  Settings":  GOSUB  7)00 
668C-  At  *  "horizontal  scale" :LB*  =  "LHz'div' :  GOSUB  7000 


8330  Hi  =  ‘Vesoiut  ioci"  :L&4  —  -  :  --.-iBUS  7^00 

8640  OS  =  "position  of  dot  ma'-Ler1  :.B*  =  "WH:"  :  GOSUB  7000 

8650  OS  =  "vertical  scale"  :lB*  =  "rc.crovolts/div":  GOSUB  7000 

8660  GOSUB  7300: AS  =  "10400  settmos:":  GOSUB  7100 

8670  PS  =  "bloct*  sire"  :  LBS  =  "tvtei":  GOSUB  7000 

868 0  PS  =  "sample  i ntervsi " : LP»  =  "mS":  G0SUB  7000 

8630  OS  =  "  Ente--  trie  vertical  range  settings":  G05UB  7100:01*  -  "Cnanr, 
el  0":03S  =  "Channel  B":lH  -  '  */—  v"  :  GOSUB  7300 
8710  ft*  =  ”  Enter  the  method  oi  coupling  -  DC/ftC"  :  G0SUB  7100:01*  =  “Ch 
arinel  0":A3'S  -  "Channel  B":LB*  =  " " :  GDSUB  7300 
8730i  AS  =  "  Enter-  the  Dias  :  -iC'SUB  7100:01*  =  "Channel  "  :  03*  =  "Chanr.e 
1  B”  : LBS  =  GOSUB  7300 

8750  AS  —  "tr-igger-  source"  :LB*  -  "  :  GOSUB  7000 

8760  ftS  =  "slope" :CHS  =  GOSUB  7400 

8770  OS  =  "trigger  coupl i ng ‘  : C-HS  =  "AC/DC":  GOSUB  7400 

8780  ftS  =  "  level  ":LBS  =  "■/■":  GOSUB  7000 

8790  OS  =  “window"  : UBS  *  GOSJB  7000 

8795  OS  =  "delay" :LBS  =  " " :  GOSUB  7000 

8800  GOSUB  7300: A*  =  "Miscellaneous  items:' :  GOSUB  7100 
8805  ftS  *  "water  pressure" : LBS  =  "osi":  GOSUB  7000 

8810  OS  =  “  Enter  pea)<-pea)<  voltages  on  transducers  1  and  3.  ":  GOSUB  71 
00:ftlS  =  "Transducer  i":A3*  =  "Transducer  £":LBS  =  "v”:  GOSUB  7300 


8830 

OS  =  “depth  of  water" :LBS  = 

"cm":  GOSUB  7000 

8830 

OS  =  "voltage  attenuation":.. 

BS  =  "  (3'0dB)  "  :  GOSUB  7000 

8840 

OS  *  "translat  ion"  :cBS  =  "  1  r. 

.  " :  GuSUB  7000 

8850 

OS  *  "direction  of  trar.siati 

on' : LBS  *  " Ri ght /Lef t " :  GOSUB 

8660 

GOSUB  7300:  HOME  :  RETURN 

8890 

REM  *******•*****■■*■**■■*■***■»« 

8695 

REM  8900  SUBROUTINE 

ASSIGNS  VOLUES  TO  ALPHAS 

8900 

REI*  ♦**•****«  *****«■**■►■**♦■•* 

8903 

ALPHAS ( 1 ) 

8904 

OLS (2)  - 

"B" 

8906 

OLS  < 3)  = 

"D" 

8908 

OLS  < 4 )  « 

"F" 

8910 

OLS  (5) 

"H" 

8913 

OLS (6)  = 

"  J" 

39 1 4 

Oi_S  (  7 )  = 

"0" 

8916 

OLS  <  8 )  * 

"C" 

8918 

OLS <91  * 

"E" 

8930 

OLS (10)  = 

"G" 

8933 

OLS (11)  = 

"I  " 

8934 

OLS (13)  = 

"K" 

8936 

OLS (13)  * 

M0'?N  " 

8938 

OLS  (14)  ■= 

"OBN" 

8930 

OLS (15)  « 

"OBN" 

8933 

OlS ( 16)  = 

"OPN" 

8934 

OLS (17)  = 

"OHN" 

8936 

Ol* (18)  « 

"OJN" 

8936 

OLS (19)  = 

"UON" 

8940 

OLS  ( £0 )  ■= 

"OCN" 

8943 

OLS (31)  « 

"OEN" 

8944 

OLS (33)  * 

"OGN" 

8946 

OLS (33)  = 

"OIN" 

8946 

OLS (34)  « 

“OKN" 

8950 

OLS (35)  = 

"N'>0" 

8953 

OLS (36)  = 

"NBO  ” 

6960 

RETURN 

8990 

REM  ***< 

6995 

REM  9000  SUBROUTINE 

DECIMAL  FORM  DF  MftChl  Nc  wJul: 

UNS  STEPPER  MOTOR 

9000 

REM  ***< 

130 


90 1 0 

On  f  H 

9030 

DATA 

9030 

DATA 

9040 

DATA 

9050 

DATA 

9990 

REM 

9995 

REM 

N  ACTS 

9999 

R£m. 

173, 83, 133, 1 £3, 3 
33,  1£8, 353,  1 73, 86 
193,  169, 43,  33,  166 
353, 303, 308,  336,  136 
308, 331 , 36 

********************* 

10000  SUBROUTINE  CONVERTS  INPUTTED  PULSE  LENGTH 

jAL  NUMBER  THE  COMPUTER  CRN  USE. 

-ft*-****#**-***-**--**'***-** 


i  N  T  l  A 


1  0000 

X  —  (310.35  — 

10  *  <13 

-  0)  ) 

11  /  £) 

10010 

Y  =  (  -  14.5  * 

X)  /  5 

1  0030 

2  =  <  -  14.5  - 

X )  /  5 

10030 

IF  V  )  Z  THEN 

G  =  Y 

10040 

IF  Z  >  Y  THEN 

G  =  Z 

.0050 

G  =  INT  <G) 

1  0060 

LET  K  =  5  *  G 

*  G  ♦  37 

*  G  ♦ 

at 

10070 

LET  K  -  K  /  3 

10080 

C  =  <310.35  - 

10  *  <12 

-  P)  ) 

(1  /  cl) 

J  0030 

D  =  <  -  14.5  ♦ 

0/5 

101 00 

E  =  <  -  14.5  - 

0/5 

10110 

IF  D  )  E  THEN 

H  =  D 

10130 

IF  E  >  D  THEN 

h  =  E 

10130 

H  =  INT  (H) 

10140 

LET  J  =  5  *  H 

*  H  ♦  £7 

*  H 

at 

10150 

LET  J  »  J  /  £ 

10160 

RET  URN 

10190 

REM  ******** ************ 

10195 

REM  10E00  SUBROUTINE 

PRINTS  OUT 

s 

CALCULATED  IN 

10000 

10300 

10310 

GOSUB  10000 

"UAl  Pu. 


10330 
.0330 
,0340 
i  0350 
:  0360 
14930 
14995 


,  3000 
15010 
'  303  0 


CHR*  (4 )  ;  "PR#  1  •' 
"Pulse  CTi  for  ";K:" 
"Pulse  off  for  “ ; J : 1 


CHR*  (4) 


PRINT 
3' R  I  NT 
PR  1  NT 
PRINT 

re Turn 

rttf  ******************** 

rem  15000  subroutine 
*J  HVErtAGc  !  RIGGER  AAM:PS — 

-  A  BV  hcul-  *  3  HOURS. 

R'-T  nOADD  IS  p  BuulEHn 
TsU-  0  AVcRAGE  triggers 


m  i  crosecorios’ 
raicrosecorios 


*#0'  :  REM  SHUT  OFF  PRIMER 


HLuOwS  UScR  - u  DEC. 
IF  YOU  DON’T  AVERAGE. 


Dt  I  r 
PROG 


HU / in, 

RAM 


HU  ■ 
hlJNS 


VARIABLE  HMD  El/UH-S  1  1-  JbLn 


1  ,  .  NT 
:  Z  >.  * 

l  *  *  = 

RE  -  u-  RN 


Hr.et  is  tne  motn  of  tie  aiso.av' 


=  0 


'  !  /  c. 


■mOhDD 
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H  ilr*l  hi-tecnr.i-jaes 

s  ar-e  oer*forwec 
ays. 

PRINT  -•...o'—'  lj: 

D I  !*:  F %■  ( 1  £>*• )  :  L*  I  -*  . -v  *  t  wfe . 

ui.v 

Li  ;  6  (  ;  .  s  Z^£*  >  .  •  \Zt*c  f  m  h-'« 

*  OF  SCAN  PG3I7I0?vra  — 
oi-bub  :  £>C5Jr<  t-.  - 

GG5U6  5*2*00  5  *t£s,i  GL~  P  *- 

-*b  :  £?  i  c-  -- 


••:•  l  a  2"  •_!/«♦.  i  CO  ^ooie  •*♦;-.  inert  stat  i  .c  ic 

.=?  ^iNvy.  M  loco  continues  for  ar»y  r.uwDer  of  err 


V  _l<6  .•  .  \  I-W  t  .  2">C.<*-  *  Z/<^'  ■ 


1  M_  4-^J  i  * 


;  0 

-Oft  JJ 

t  FI«s;  < 

i  c 

GOSlici  3 

000 

i  j 

Lis.  =  (  h  / 

cv 

4  *  *  CHR*  ( Co > 

1W 

PRINT 

D*; “PR* 7" 

i  tc 

PRINT 

D*  ;  "  Inm  7" 

I5id 

PRINT 

■'  SO  1  “ 

^  GO 

PRINT 

“  RA  M 

l&d 

PR  1  NT 

“AT  ■•';/*: 

1  65 

X  =  1 

190 

PRINT 

••ADA-  :  Z*: 

COO 

GET  &* 

r,CE^o  “ 


£10  IF  A3C  (b*i  =  13  Then  GOTO  £50 
£3*  IF  B*  ■  THEN  GOTO  £40 

£30  C*  *  C*  ♦  E<* 

£40  GOTO  £00 

£50  PRINT  C* : F*  < * )  =  C* 

£51  X  =  X  ♦  1 
£54  C*  ~  “  " 

£57  IF  X  *  139  THEN  GOTO  £70 
£60  GOTO  190 
£70  PRINT  "UT" 

£00  PRINT  D* ;  "F'AttO" 

£90  PRINT  D*;"IN«tO“ 

£91  REM  reset  10400  after-  cat  a  transfer. 

£3£  A  =  PEEK  <  -  16£9 1 ) :  FOR  ]  =  i  TO  1 OOO  STEP 
£33  A  -  PEEK  <  -  16£9£>:  FOR  I  =  1  TO  10000  STEP 


imE/ 

:  NE  .v  T 


£35 

FLASH 

£36 

PRINT  "WORKING 

ON 

SUBROUTINE  1 OOO " 

£97 

NORMAL 

300 

GOSU6  lOOO 

350 

FLASH 

35 1 

PRINT  "WORKING 

ON 

SUEiROUTINE  £000" 

NORMAL 

400 

GOSUE<  £000 

500 

NEXT  JJ 

510 

FLASH  :  PRINT  " 

SAVING  DATA  TO  DISK 

550 

GOSUEi  6000 

600 

END 

1000 

FOR  X  =  1  TO  1 

38 

ST  EP  1 

1005 

J  =  O 

NORMAL 


lOlO  IF  X  > 


6  AND  X  < 


63  THEN  J  =  X 


'll  flj  ”1  'XI  I'iJ  rtj  Mj  alj  .;u  "J  Mj  |1J  Mj 


A  \ 


IF  J  —  0  THEN  GOTO  1  1l5 
L  -  1 

K  =  ( 6  *  J  )  - 

H*  =  MID*.  (FS(X;  .  3,  13) 

G  (  K  )  —  VAl  ( H*»  ) 

M  =  14 

FOR  i-  =  3  TO  B  BfEP  1 
K=  (6  *  J  >  -  (  6  -  L  > 

H*  =  MID*  ( F*>  t  a  >  .  M,  1  1  > 

G(K>  =  VA  l  (H*) 

M  =  M  ♦  1  1 

NEXT  l 
NEXT  X 
RETURN 

REM  STATISTICAL  CALCULATIONS 
REM  IF  CONTINUOUS  iNPl)T  OF  P 
Cl  =  16. 7 j  REM  Cl  AND  03  ARE 
NT  ALL V  DETERMINED 
Cl  —  j. Sll 

FOR  K  *  P  TO  0  STEP  1 
I  (K)  =  (Cl  *  G(K)  +  C3)  3 

NEXT  K 


OF  P'  AnO  G1  NlCESSARY  GOBuB  5000  hERE 

ARE  CALIBRATION  FACTORS  THAT  WERE  EXPERImE 


FOR  K  =  P  1 
XX  =  <1  /  L ) 

AA  *  AA  *  X> 
SB  *  BB  ♦  X) 
NEXT  K 
II < JJ)  *  AA 
FF(JJ)  =  BB 


C'  STEP 
►  (  1  i  K  - 


:  DI)  =  (?:  EE  =  0 
=  P  TO  0  STEP 


Li7e. 

XX  =  (1 

/ 

3) 

« 

i  1  ( K  *  1  )  * 

I  IK)  ) 

3  l  80 

CC  =  CC 

+ 

XX 

* 

( K  —  FF  *  J  J )  ) 

2 

3  a  30 

DD  =  DD 

■+ 

XX 

* 

< K  -  FF  ( JJ)  ) 

5 

3301? 

EE  =  EE 

XX 

* 

(K  -  FF ( JJ> ) 

4 

3310 

NEXT  K 

3330 

VAR ( JJ) 

= 

CC 

/ 

I I ( JJ) 

3330 

RMS ( J J  > 

= 

SQR 

(VAC 

3340 

SKEW ( JJ) 

= 

DD 

/ 

< 1 I < JJ)  *  (RMS ( JJ) 

2) 

3350 

KURT ( JJ ) 

= 

EE 

/ 

i  I  I  ( JJ)  *  (RMS  •;  JJ  < 

4  > 

3360 

PRINT 

D* :  "PR# 1  " 

3365 

PRINT 

"SCAN  POS=  " 

JJ 

3370 

PRINT 

"P  CURSOR*  " 

p 

3360 

PRINT 

"0  CURSOR* 

C' 

3385 

REM  Xr-  IS  THE  VUcTAGE  ATTEfnuA 
COMPUTE  INTENSIFY 

3330 

PRINT 

" INTENSITY* 

;  I  I  i  J  J  )  /  ' 

3300 

PRINT 

"AVE  FRED*  " 

FF  ( JJ> 

3310 

PRINT 

"RMS  FRED*  " 

RMS ( JJ) 

3330 

PRINT 

"VARIANCE*  " 

VAR ( JJ> 

l330 

P'  A  I  N  T 

"SKEWNESS* 

SKEW ( JJ ) 

3340 

PRINT 

"KURTOSIS*  ■ 

KURT  < JJ) 
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£250  PRINT  D*;"Pft*0' 

£260  RETURN 
3000  Ds  =  CHRS  <4) 

3010  i*  =  CHRl  (£t,) 

30£0  P'RINT  D"*:"P'R87" 

30.30  P'S  IN  T  Hi  ;"IN#7' 

3040  P'RINT  '  .-jC  1  " 

3050  P'RINT  ’  wp" 

3060  £l»  =  PPP "  ♦  NPM*  ♦  "  =  " 

30  70  E3*  =  "=PPG" 

3060  PS  =  STRt  (JJ) 

303f0  £*  =  El*  ■*  PS  +  E£* 

3100  P'RINT  "UT'";ZS;E* 

3110  FClR  I  =  1  TO  15000  STEP’  1:  NExT  I 

3115  PRINT  "UT" 

3116  PRINT  D*;"PR#0" 

3117  PRINT  D*;"IN#0" 

31£0'  P  =  PEEK  (  -  16331; 

3130  FOR  I  =  1  TO  1000  STEP  1:  NEXT  I 
3140  P  =  PEEK  (  -  16333) 

3:50  FOR  I  =  1  TO  10000  STEP  1:  NExT  I 

3160  RETURN 

4000  REM  INITIRl  DPTP 

4010  SuSUR  6300:  REM  GET  T  I Tl_E  PnD  KEYSTROKES  FOR  TJ 

4030  PRINT  :  P'RINT  "Wnat  is  the  first  scan  DOSit  lcr." :  :  InP_"  ;  <S1 

4030  PRINT  :  INP'UT  "Nhat  is  the  last  scan  position  ":lPS~ 

4040  INP'UT  "Enter  voltage  attenuation  factor  in  dR  ";XF 

4050  PRINT  CHRS  <4>;"PR#1" 

4055  INPUT  "Enter  the  angular  setting  of  RPLMIKE  hut  :j.,s  :  3 

4060  P'RINT  "Data  from  scan  positions  FIRST;"  to  :lPST 

4065  IF  LPST  -  FIRST  >  50  THEN  PRINT  "ERROR1  •  This  program  only  ac 

cept  50  scan  positions.  To  input  more  cnanpe  the  a i mens i pts  of  the 
variaoles  on  line  6  and  this  line,  4065.  5TQP' 

4070  P'RINT  "*  NOTE:  Attenuation  factor  of  "  ;  XF  ;  "  d  is  usee  t  .  uchout 
calculations. " 

407!  PRINT  "  "  ;  X  F  ; "  dE1  translates  to  a  decimal  factor  ::  GOSuR  90 

00:  PRINT  XF:  REM  3000  CONVERTS  XF  TO  DECIMP;.  VALUE 
4075  P'RINT  ”  Plso,  tne  angular  setting  of  RPLMir.E  is  :'3:  oeg-'ees. 

4077  PRINT  :  PRINT 

4080  PRINT  CHRS  <4>;”PR#0" 

4030  RETURN 

5000  REM  THIS  ROUTINE  MINDS  ITS  P’s  PND  O'  s 

5005  REM  P  PND  Q  MPRK  THE  LEFT  PND  RIGHT  ROUNDS  I E6  TO  ;.VE;,c,L':  £_C-  L. 
RPY 

5010  PRINT  "SCPN  POSITION  =■  ";JJ 

5040  PRINT  :  P'RINT  "Enter  P  "  :  :  INPUT  P». 

50S0  P'RINT  "Enter  0  INP'UT  C'S 

5060  IF  LEN  (PS)  =  0  THEN  5040:  REM  DONT  CHPNGr  F 

5070  P  «=  VPL  IP'S) 

5030  IF  LEN  (OS)  =0  THEN  5040 

5030  G  *  VPL  (C*i 

5100  IF  0  <  =  P'  THEN  P'RINT  "ERROR  '  '  f'  neecs  *  ma.-M  r  it  —  •  ■  *  oc1 

rtion  of  the  domain  and  0  the  right. GOTO  5000 
5110  RETURN 

6000  REM  SPVE  DPTP  TO  D 1 3k 

6005  IF  (  PSC  (VIS)  >  47)  PnD  i  HuC  i  T  IS  i  .  v' 

"Since  1st  character  of  title  is  a  P.  :  -win. 

60  10  F  ‘  I  7  LES  —  TIT  l£*  ♦  .  I  NT  EL N 3  i  !  i  1 


.  s  : 


6030  PRINT  D*; "WRITE  •  ;  F';  I u£l 

6035  REM  LAST  =  lAS-'  SCAN  POS 1  T  i  UN* 

6037  REM  rIRS:  -  FIRST  5CR  w  .-‘uSJ  i  l)n 

6040  PRINT  <_AS  .  -  ■  I  RS  I  »  1 

6050  FOR  JJ  =  -IRS"  TO  ..AST  =ltp  ; 

6060  PRi.'vT  Jj 

6070  PRINT  lluji 

6060  Tt>  i  JJ 

6030  PRINT  D*  :  "LLUoc: 

o<-3p  L»*  ;  _ULr.  ;F  I  i  :=-3* 

b  1 00  FTITLE*  -  I  TlE*  ".AvEftWoE  FRED" 

o.l.E'0  PRINT  D*;"GPEn  ;  FT  1  T  Lc* 

6130  PRINT  D*; "WRITE  "jFTITlE* 

6140  print  (last  -  fi rstj  ♦  i 

£150  FOR  JJ  •=  FIRST  T0  i_AS~  STEP  . 

£160  PRINT  JJ 
6170.  PRINT  FFlJJ> 
bi60  NEXT  JJ 
6130  PRINT  D*!"CLGb£ 

6i  35  PRINT  D* ;  "LOCr  7l.E* 

6300  FTITlE*  =  TlVuE*  -  .  HMS  FRED" 

6l30  PRINT  Q*  :  uPcN  jFTITle* 

6330  PRINT  0*  ;  "wRI E  FTITlE* 

£340  PRINT  ;lAS.  F  *  i 

6350  FOR  jj  =  FIRS’,  to  last  STEP  1 

6c.  60  -Ki.vr  jj 

6370  PRINT  RMSiJJ) 

=360  Nfc  A  JJ 

6330'  PRINT  D*;  "CLOSE" 

o-3p  PRLNI  09  ;  "  LOCr,  FTITlE* 

6300  FTITLE*  =  TITLE*  *■  ".VARIANCE 
6330  PRINT  D*;  OPEN  ;FTITLE* 

6330  PRINT  D*;'wRI~E  ":FTITl£* 

£340.  PRINT  (LAST  -  f  ;  Si  ,  )  1 

6350  FOR  JJ  *  FIRS’  TO  ..AST  STEP  I 
6360  PRINT  JJ 
6370  PRINT  VAfiUJ) 

6360  NEXT  JJ 

6330  PRINT  D»: "CLOSE 

6335  PRiNT  D*;"LOCl,  "  ;FTITLE* 

6400  ft  i  T  i_E*  =  TITLE*  *  ".SKEWNESS" 


o4,£0 

PRINT  D*  ;  " OC'EN 

" I TL6* 

6s30 

PRINT  09; "WRITE 

*'  ;  PT  I  Ti_Ef» 

t >4*0 

PRINT  (LAST  -  f 

I kS7  >  ■+■  1 

6451? 

FOR  JJ  =  FJRS’: 

TO  last  STEP  i 

6460 

PRINT  Jj 

647i? 

PRINT  SKEW(JJ) 

r»40?i 

n£  x  f  j  J 

6490 

PRINT  D*; "CLOSE" 

b*+ 

PR  I  >iT  D*  ;  "  lUlk 

;  F  i  i  * 

6500 

F  T I TLE*  =  TITLE* 

. KURT0S1S 

6  5  JO 

PRINT  D*;"OPE. 

iFTITlE* 

6530 

PRINT  D*; "WR j 1 E 

"  :FTITlE* 

65-40 

PRINT  (LAS!  -  F 

l  AST  >  *  1 

6550 

FOR  JJ  =  FIRS'. 

TO  LAST  STEP  1 

6560 

PRINT  jj 

6570 

PR i NT  KJA'  IJJ- 

bliUv 

I*-  X  i  J  j 

65  s0 

PRINT  D*;'C_CoE 

C/i<  -ft- 

PR : N i  D* :  _OCr 

"  jFTITLE* 

6600 

RETURN 

137 


ut:  ;  rt  Trri-t:  *04  Dw  lh 

6305  =  " " 

a3iC  PR  INI  :  INPUT  "Enter  the  title  or  the  fi  ;e  to  oe  read  "  : 

S315  IF  LEN  (TITLE*)  >  5  TrtEN  PRINT  'ERROR-  Tne  maximum  length  or  a 

oats  laDel  is  5  cnaractet’5.  " :  GOTO  6305 

8317  IF  lEN  (-TITLE*)  (  1  THEN  PRINT  "ERROP'  1  You  must  enter  sc-rnethinri 
.  "  :  GOTO  8305 

8380  FOR  J  =  1  TO  lEN  t  TITLE*) 

S33&  CHAR*  =  MID*  (TITLE*.  J,  1  >  :  HEM  GET  A  CHARACTER  u- r  0  F  STRING 
83<*0  IF  (  ASC  (CHAR*)  )  36)  AND  (  ASC  (CHAR*)  <  183)  TrtEN  CHAR*  =  CHR* 

(  ASL  (ChAR*)  -  38)  :  REM  MAKES  LOWER  CASE  INTO  UPPER  CASE 
835®  IF  (  ASC  (CHAR*)  >  64)  AND  (  ASC  (CHAR*)  (  31)  THEN  NCHAR*  =  AuPhAj 

(  ASC  (CHAR*)  -  64) :  GOTO  838®:  REm  IF  LETTER  THEN  DECODE 
636®  IF  (  ASC  (CHAR*)  )  *.7)  AND  (  ASC  (ChAR*>  <  58)  THEN  NCHAR*  =  ChAR*: 

GOTO  836®:  REM  A  NUMBER  HAS  SEEN  CHOOSEN 
837®  PRINT  :  PRINT  "ERROR1 •  You  can  only  enter  a 1 phariumer I cchar act ers. " 
:  GOTO  63®5 

638®  NAM*  =  NAM*  ♦  NCHAR* 

639®  NEXT  J 

8393  PRINT  CHR*  (4) ; "PR*1 " :  PRINT  "Data  file  name. . . " s T I TLE* 

8394  PRINT  CHR*  ( 4 )  ;  "Pfi*=® " 

6395  RETURN 

89®®  REM  ASSIGN  MACROS  TO  LETTERS 
8908  ALPHA*(1)  =  "i?" 

6904  Al* (8)  =  "6" 

8906  AL* (3)  =  "D" 

8906  AL*  ( 4 )  =  "F" 

8910  AL* (5)  =  "H" 

8318  Al* (6)  =  "J" 

8914  AL* ( 7 )  =  "A" 

39 l 6  AL* (6)  =  "C" 

8916  Al* (9)  =  "E" 

6960  AL* (10)  *  "6' 

8388  AL*  (  i  1  )  =  "  I  " 

6384  At.*  (  18)  =  "K" 

6986  Al*  (13)  =  "OC»N" 

8388  AL* ( 14t  =  "OBN' 

6930  AL* (15)  *  "ODN" 

8338  Al* ( 1 6  >  =  "OFN ' 

8934  AL* (17)  =  "GHN " 

8336  AL *(18)  =  "OJN" 

8338  AL*  < 1 3 1  =  ”OAN" 

834®  AL* ( 80 )  =  "OCN" 

8348  AL* (81)  =  ”OEN' 

8344  AL  *  ( 88 )  =  "OGN 

8946  AL* (83)  =  'OIN' 

6346  Al*(84j  =  'OKN' 

8350  Al*  ( 85 )  *  "NC«0' 

8358  AL*  <  86 )  =  "NBO" 

836®  RETURN 

9000  REM  db  TO  DEC! Mho  CONVERTER 

3010  XF  =  -  1  *  XF  /  80:  REM  BECAUSE  VOi-TAGE  IS  THE  SORT  OF  POWER 

9080  XF  =10  XF  :  REM  CONVERTS 
9030  RETURN 
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{note:  no  decimals  in  the  angles} 

{updated  25Mar92} 
program  Macintosh_analyzer; 

uses 
dos,  crt; 

const 

maxpts  =  40;  {equal  to  max  angles  any  one  data  run  will 
have} 

{unlikely  to  be  larger 

than  20} 

offset  =  3.823;  {updated  18Mar92  to  25Feb92  calibrations} 
vdiv  =  19.8162; 

interval  =  0.050;  {seconds  per  point} 
noisefactor  =  0.00;  {noise  contribution  per  point} 
path  =  'a:\'; 
type 

v_ary  =  array  [1  ..2000]  of  real; 
str80  =  string  [80]; 
str40  =  string[40]; 
str20  =  string[20]; 
str8  =  string[8]; 

realArry  =  array  [1.. maxpts]  of  real; 
strAry  =  array  [1..  10]  of  str8; 

var 
f  :  text; 
volt  :  v_ary; 

NumFiles,  i,  runs,  last_pt,  p,  q,  run_no,  dBs  :  integer; 
x,  time,  angle,  current_angle,  actual_angle,  dBFactor  :  real; 

ExperNum,  dataname,  possibly  :  str8; 

PtVal,  Xavg,  Variance,  Rms,  Skev/ness,  Kurtosis,  intensity2, 
time_irtens2f  thirdmom,  fourthmom  :  realArry; 
nextFile  :  SearchRec; 
fil  :  str80; 

fileName,  ExperAry  :  strAry; 


procedure  getFilelnfo; 
var 

i  :  integer; 


1  40 


begin 

write('Number  of  files  to  read?  '); 
readln(NumFiles); 
for  i  :=  1  to  NumFiles  do 
begin 

write(’Enter  name  of  file  ',i,'  ’); 

readln(fileName[i]); 

write(’Enter  Experiment#  of  file  ',i,'  '); 

readln(ExperAry[i]); 

end; 

end; 


procedure  nextAngle(var  A  :  integer); 
var 

si  :  str80; 
error  :  integer; 


begin 

FindNext(nextfile); 
if  doserror  =  0  then 
begin 

si  :=  copy(nextfile.name,6,3);{ strip  out  angle) 
val(sl.  A,  error); 
end 

else  {errors  mean  last  file  read) 
a  :=  -99; 

end; 

function  convert_angle(theta  ;  integer)  :  integer; 
begin 

theta  :=  theta  -  284; 
if  theta  <  -180  then 
theta  :=  theta  +360; 

convert_angle  :=  theta; {sets  range  at  -180  to  180} 
end; 

procedure  getFirstFile(var  a  :  integer); 
var 


s2  :  str80; 
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sl  :  str20; 
error  :  integer; 

begin 

s2  :=  copy(dataname,  1,4); {strip  4  characters} 
s2  :=  concat(s2,ExperNum); 

s2  :=  concat(path,  dataname.XJ'.ExperNum,'^',  s2,'*.asc’); 
repeat 

findFirst(s2,  anyfile,  nextfile); 
if  doserror  o  0  then 
begin 

writeln('Check  to  see  if  ',path,'  has  correct 

data  disk.'); 

writeln(T'm  looking  for  file  ’,s2); 
repeat  until  readkey  =  #13; 
end; 

until  doserror  =  0; 

sl  :=  copy(nextfile.name,  6,3); 

val(  sl,  a,  error); 

end; 

procedure  alternateparser(al  :  str80;  var  x  :  real); 
var 

error  :  integer; 


begin 

val(al,  x  ,  error); 

end; 

procedure  readwave(filenam  :  str80); 

var 
f  :  text; 

pt  :  integer; 
temp_str  :  str80; 

begin 

pt  :=  0; 

assign(f,  filenam); 
reset(f); 

readln(f,  temp_str);  { first  2  lines  are  useless  to  me} 
readln(f,  temp_str); 


while  not  (eof(f)  or  (pt  >  1999))  do 
begin 
pt  :=  pt  +1; 
readln(f,  temp_str); 

temp_str  :=  copy(temp_str,  5,  100);  (strip  first  4 
characters  from 

temp_str  --  ie.  "  "  } 

altemateparser(  temp_str  ,  volt[pt]); 
end; 
close(f); 
last_pt  :=  pt; 
end; 

procedure  get_pq; 


filenam  :  str40; 
f  :  text; 

angle,  file_p,  file_q  :  real; 

begin 

filenam  := 

concat(’b:\dataname,,_\ExperNum,\\dataname,’_,,ExperNum,’.pq 

m*); 

assign(f,  filenam); 

reset(f); 

repeat 

readln(f,angle,  file_p,  file_q); 
writeln('angle  =  ’.angle:  12:5,’  filep  =  \file_p:12:5); 
until  (trunc(angle)  =  run_no)  or  (eof(f));{  again  this  is  an 
approximation} 

{that  works  only  for  half  and  integer  angles} 
close(f); 

if  trunc(angle)  <>  run_no  then 
begin 

writelnCError!!  could  not  find  angle  \run_no,’  in 

.pqm  file'); 

q  :=  -1; 
end 

else 

begin 

p  :=  round(file_p); 
q  :=  round(file_q); 


function  dB_to_decimal(dB:integer):real; 
var 

exponent  :  real; 

begin 

exponent  ;=  (-1)  *  dB  /  20  *  ln(10); 
dB_to_decimal  :=  exp(exponent); 
end; 

procedure  getdBatten(var  attenuation:  integer); 
var 

filenam  :  str40; 
f  :  text; 
angle  :  real; 

begin 

filenam  := 

concat('b:,,dataname,'_',ExperNum,\,dataname,’_,,ExperNum,'.db 

x'); 

assignff,  filenam); 
reset(f); 

readln(f,  attenuation); 
close(f); 

dBFactor  :=  db_to_decimal(attenuation); 

end; 


procedure  calibrate; 


i  :  integer; 

begin 

writelnfstarting  calibration'); 
for  i  :=  p  to  q+1  do 
begin 

volt[i]  :=  vdiv*volt[i]+offset; 
if  volt[i]  <  0  then 
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begin 

writeln(’underflow  occured  at  pt  \i); 
end; 

volt[i]  :=  volt[i]*volt[i]; 
end; 

end; 

procedure  integrate l(var  intsum,  averval  :  real); 
var 

i  :  integer; 

xx,  firstmom  :  real; 


begin 

writeln('starting  integration  1  ’); 
intsum  :=  0; 
firstmom  :=  0; 
for  i  :=  p  to  q  do 
begin 

xx  :=  1/2  *(volt[i]  +volt[i+l]); 
intsum  :=  intsum  +xx; 
firstmom  :=  firstmom  +  xx*i 
end; 

averval  :=  firstmom  /  intsum 
end; 

procedure  integrate2(intensity,  avgval  :  real;var  Variance, 
Rms,  Sk, 

third,  fourth,  Kurt  :  real); 
var 

i  :  integer; 

xx,  cc,  dd,  ee  :  real; 

begin 

writeln(’starting  integration2'); 
cc  :=  0; 
dd  ;=  0; 
ee  :=  0; 

for  i  :=  p  to  q  do 
begin 


xx  :=  1/2  *(volt[i]  +volt[i+l]); 
cc  :=  cc  +  xx  *  sqr(i  -  avgval); 
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dd  :=  dd  +  xx  *  sqr(i  -  avgval)  *  (i  -  avgval); 
ee  :=  ee  +  xx  *  sqr  (  sqr  (i-  avgval)); 
end; 

variance  :=  cc  /  intensity; 

rms  :=  sqrt( variance); 

third  :=  dd/intensity; 

fourth  :=  ee  /  intensity; 

sk  :=  dd  /  (intensity  *  rms  *  rms  *  rms); 

kurt  :=  ee  /  (intensity  *  sqr  (sqr  (rms))); 

end; 

procedure  noise_reduc(var  noise  :  real); 
begin 

noise  ;=  (q  -  p)  *  noisefactor; 
end; 

procedure  savedata(lastposit  :  integer); 
var 

f  :  text; 

i  :  integer; 
filenam  ;  str40; 

begin 

{save  average  frequency} 
filenam 

concat('b:\dataname,'__\ExperNum,\\dataname,’_\ExperNum,’.av 

g’); 

assign(f,  filenam); 
rewrite(f); 

for  i  :=  1  to  lastposit  do 
writeln(f,  ptvalfi],  *  ’,Xavg[i]); 
close(f); 

{save  variance) 
filenam  := 

concat(’b:\dataname,,__\ExperNum,\\dataname,'_\ExperNum,'.va 

r’); 

assign(f,  filenam); 
rewrite(f); 

for  i  :=  1  to  lastposit  do 
writeln(f,  ptvalfi],'  Variance[i]); 
close(f); 


{save  rms  frequency} 
filenam  := 

concat(,b:,,dataname,,_',ExperNum,\\dataname,’_',ExperNum,,.r 

ms'); 

assign(f,  filenam); 
rewrite(f); 

for  i  :=  1  to  lastposit  do 
writeln(f,  ptval[ij,’  Rms[i]); 
close(f); 

{save  skewness} 
filenam  := 

concat(’b:,,dataname,'_,,ExperNum,'\\dataname,,_,,ExperNum,,.sk 

w'); 

assign(f,  filenam); 
rewrite(f); 

for  i  :=  1  to  lastposit  do 
writelnff,  ptvalfi],'  Skewness[i}); 
close(f); 

{save  kurtosis} 
filenam  := 

concat('b:\dataname,,_\ExperNum,\\dataname,,_,,ExperNum,’.kr 

f); 

assignff,  filenam); 
rewrite(f); 

for  i  :=  1  to  lastposit  do 
writeln(f,  ptval[i],'  Kurtosis[i]); 
close(f); 

{save  intensity} 
filenam  := 

concat(,b:',dataname,'_',ExperNum,\\dataname,’_,,ExperNum,,.i2’ 

); 

assign(f,  filenam); 
rewrite(f); 

for  i  :=  1  to  lastposit  do 
writeln(f,  ptval[i],'  intensity2[i]); 
close(f); 

{save  third  moment} 
filenam  := 

concat('b:\dataname,'_\ExperNum,\\dataname,'_',ExperNum,’.3r 

d’); 

assign(f,  filenam); 
rewrite(f); 

for  i  :=  1  to  lastposit  do 
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writeln(f,  ptvalfi],’  thirdmom[i]); 
close(f); 

{save  fourth  moment} 
filenam  := 

concat(’b:,,dataname,'_\ExperNum,'Y,dataname,'_\ExperNum,'.4t 

h’); 

assign(f,  filenam); 
rewrite(f); 

for  i  :=  1  to  lastposit  do 

writeln(f,  ptval[i],‘  fourthmom[i]); 

close(f); 

end; 

procedure  define_intensity(integral  ;  real;  var  intensity2,  time2 
:  real); 
var 

xx,  timeintvl  :  real; 

begin  {two  methods  are  used  to  calculate  intensity  because 
I'm  not  sure 

where  noise  is  } 
noise_reduc(xx); 

intensity2  :=  (integral  /  (dBFactor  *  dBFactor))-  xx; 
timeintvl  :=  (q  -  p)  *  interval;  {converts  points  to  time 
domain} 

time2  intensity2  *  timeintvl; 
end; 


(♦main  program*) 
begin 

writeln(,Don"t  forget  to  update  noiselevel  in  the  constant 

section!!!'); 

write('Stop???'); 

readln(possibly); 

if  upcase(possibly[l])  =  'Y'  then  exit; 
getFilelnfo; 

for  runs  ;=  1  to  NumFiles  do 
begin 

dataname  :=  fileName(runs]; 

ExperNum  :=  Exper  Ary  {runs}; 
getFirstFile(run_no) ; 
getDBAtten(dBs); 
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if  dBs  =  -100  then  (error  occurred) 
run_.no  :=  -100;  (and  stop) 
i  :=  0; 

while  run_no  >=  0  do 
begin 

fil  :=  concaKpath.dataname.’C^.ExperNum/jY,  nextfile.name); 
readwave(fil); 
get_pq; 
calibrate; 
i  :=  i  +  1; 

integratel(x,  Xavg[i]); 

integrate2(x  ,  Xavg[i],  Variancefi].  Rms[i],  Skewness(i), 
thirdmom[i],  fourth  mom  [i],  Kurtosis(iJ); 

define_intensity(x,intensity2[i],time_intens2[i]); 
ptval[i]  :=  convert_angle(run_no); 
writeln('At  angular  position  *,ptval[i]); 
writeln(’The  runNumber  is  ’,run_no); 
writeln('The  result  of  integration  is  \x:12:5); 
writeln('Intensity  information:’); 
writelnO  Integral  *  dB’); 

writelnO  Intensity  =  \intensity2[i]:12:5); 

writeln('  w/  time  interval  =  \time_intens2[i]:12:5); 
writeln('The  average  frequency  =  ’,Xavg[i]:12:5); 
writeln('The  variance  =  ’,Variance[i]:12:5); 

writeln(’The  Rms  frequency  =  ’,Rms[i]:12;5); 
writeln('The  Skewness  =  \Skewness[i]:12:5); 

writeln('The  Kurtosis  =  \Kurtosis[i]:12:5); 

writeln('The  Third  moment  =  \thirdmom[i]:12:5); 

writeln('The  Fourth  moment  =  \fourthmom[i]:12:5); 

writeln; 

writeln('Note:  intensity  and  average  frequency  values  are 
subject  to'); 

writeln('other  calibrations  that  have  not  been  included  in  the 
current'); 

writeln('calculations.’); 

writeln; 

NextAngle(run_no); 

end; 

if  run_no  <>  -100  then  (might  have  problem  with  dbx  file) 
savedata(i); 
end;  (for  runs) 


